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CHAPTER I 


INTRODUCTION 

In his Treatise on the Mathematical Theory of Elasticity, 

Love* observed that "When the general equations had been obtained, 
all questions of the small strain of elastic bodies were reduced to 
a matter of mathematical calculation." To this day, that 'matter of 
mathematical calculation' figures prominently in applied mechanics. 

The early mecanicians realized that the general equations of 
elasticity were too difficult to solve except in a few special cases, 
so a large part of their effort was focused on methods for finding 
approximate solutions to problems of technological interest. Some of 
the techniques they used in deriving approximate theories for rods, 
plates, and shells are, in fundamental ways, very similar to the 
finite element technique. 

Today it is well understood that the classical theories of 
rods, plates, and shells may all be systematically derived from 
elasticity theory by introduction of approximations for the displace- 
ment to the principle of virtual work. Kirchhoff is the first person 
mentioned by Love [1] as having used this methodology,** and in using 
it Kirchhoff managed to give a clear interpretation of the boundary 

Love [1], p. 2. 

.C A 

to derive a theory of elastic rods; later, elastic plates. 
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conditions in the plate theory with which his name Is now associated. 

In two respects Kirchhoff's methodology is the same as the finite 
element methodology. First of all, he made kinematic approximations, 
and secondly, he used an energy principle to maintain consistency 
between his generalized stresses and strains, and to arrive at the 
correct boundary conditions. The principal difference between Kirch- 
hoff's methodology and the finite element methodology lies in the 
degree to which the kinematic field is approximated. Because of the 
similarities in the construction of the classical rod, plate, and 
shell theories to the construction of finite element equations, the 
successes and failures of the classical structural theories reflect, 
at least qualitatively, upon the performance of the finite element 
method. 

No special theory in the realm of solid mechanics has enjoyed 
greater success than that of elastic beams, for there the general 
equations of elasticity are effectively replaced by a single ordinary 
differential equation. The theory is not only reasonably accurate, 
but extremely easy to understand because of its displacement based 
derivation. The classical plate and shell theories provide equations 
less easy to understand and less easy to solve than the beam equations, 
but still regarded as simpler than the general equations of elasticity. 

A major failing of the classical theories of beams, plates, 
and shells is their Inability to account for the effects of 'trans- 
verse shear stress'; that is, the shear stress acting on plane 
sections through the thickness of the structure. As a direct con- 
quence, those theories always give a higher estimate of the stiffness 
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of a structure than does the general theory. Secondly, the twisting 
moment and shear force are coupled on the edge of such a plate or 
shell. In spite of these shortcomings, it was not until after 
Reissner's [2] investigation into the effect of shear stress on the 
bending of plates that satisfactory alternatives to the classical 
theories became available. But Reissner's methodology has had a 
greater impact on the methods used in applied mechanics than did his 
plate theory of itself. In its derivation his theory is distinguished 
from the classical theories by the fact that both assumed stresses and 
displacements are used. Since that time the use of assumed stresses 
in the derivation of plate and shell theories has become common. 

It is not surprising that the finite element method has evolved 
along similar lines. The mot i ve--f i ndi ng approximate solutions for 
problems of technological interest--was the same for the early finite 
element researchers as it had been for the early plate and shell 
theorists. The finite element method in which one introduces kinematic 
approximations to the virtual work principle is the direct counterpart 
of Kirchhoff's rod and plate theories. The same types of advantages 
and defects are inherent. 

The principal advantage of the displacement based finite 
element methods is their conceptual simplicity. For application to 
beams, the simplicity rivals the simplicity of the beam theory itself. 
In the cases of plates and shells though, it proves difficult to con- 
struct 'compatible' shape functions for the displacement. Finite 
elements for thin plates based on kinematic approximations sometimes 

overestimate the stiffness of the plate so badly that they are 
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described as 'locking.' As a means of avoiding locking, and just for 
simpler construction of shape functions, some researchers have pre- 
sented 'incompatible' plate bending elements, elements which do not 
satisfy slope continuity at element interfaces. A second problem with 
displacement based finite element methods (in general) is their ability 
only to satisfy traction boundary conditions In an average sense. One 
of the principal advantages of the finite element method over the 
method of finite differences is its ability to satisfy higher order 
boundary conditions accurately; but this potential is not fully 
realized in a purely kinematic formulation. 

It was Pian's [3] investigation into the derivation of element 
stiffness matrices that brought widespread attention to the potential 
advantages of introducing the stress as an independent variable. By 
his formulation, which was based on the complementary energy principle 
of linear elasticity, an energy-consistent alternative to incompatible 
elements was made available. Also, as was the case in Reissner's 
plate theory, the stress formulation made possible considerably more 
accurate satisfaction of traction type boundary conditions. Finally, 
Plan observed a marked acceleration in the convergence of the 
components of the stiffness matrix when the stress method was used. 
Since that time, the study of fini te e lement methods related to 
Pian's (which have come to be known as 'hybrid stress methods') has 
produced a number of special methods which may be applied where con- 
ventional displacement based finite elements fail. 

The failure of the conventional finite element method seems 
always to be attributable to the presence of kinematic constraints 
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which cannot be satisfied by approximated displacement fields. In the 
case of plates and shells, that constraint is interelement slope 
continuity. When the algorithm is based on an energy principle, it is 
always possible to 'relax' the constraint by the introduction of a 
Lagrange multiplier. That multiplier always will be a generalized 
stress, the 'energy-conjugate' of the kinematic constraint. 

One particular class of problems in which the conventional 
displacement-based finite element method fails is composed of problems 
Involving incompressible or nearly incompressible bodies. The 
constitutive equation for such bodies is nearly or precisely singular 
for the mode of dil Station. The shape functions for the displacement 
used in the conventional finite element method are incapable of pro- 
ducing any motion other than pure shearing which does not contain 
(loosely speaking) 'excessive' dilitation. As a consequence, the con- 
ventional finite element method drastically overestimates the ' 
resistance to deformation of nearly and precisely incompressible 
bodies. In a key paper by Herrmann [1*], it was shown that the dif- 
ficulty could be avoided if only the mean stress were introduced as 
an independent variable. In this case the kinematic constraint was 
incompressibility and the energy-conjugate stress was the pressure. 

His assumed pressure formulation of finite elements for such bodies 
is now o standard practice. 

Problems involving finite deformations of strain-softening 
bodies resemble problems involving nearly incompressible bodies In the 
sense that the body's shear compliance is much greater than its bulk 
compliance. For the most part, finite element analyses of such 
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deformations have been accomplished only at considerable expense, even 
when the pressure is introduced as an independent variable. No finite 
deformation counterpart to the complementary energy principle of linear 
elasticity was known, so no pure stress or hybrid stress finite element 
algorithm was found. 

The door to stress based finite element analysis of f TnTt § 
deformation problems was opened in 1972 by Fraeijs de Veubeke [5] with 
his presentation of a complementary energy principle for finite defor- 
mation elasticity.* The stationary conditions of this principle are 
both the equations of compatibility and angular momentum balance. To 
date, variants of the principle have been used by de Veubeke and 
Millard [6], Sander and Carnoy [/l, Koiter [8], Wunderlich and 
Obrecht [9], Murakawa [10], Murakawa and Atluri [11], [12], Murakawa 
et al. [13], and Atluri and Murakawa [1A], in problems ranging from 
elastic membrane theory to beam, plate, and shell theories. 

A considerable generalization of de Veubeke's principle was 
given by Atluri [15]. His formulation of the complementary principle 
for stress rates and spin opened the way for the current work, that 
of developing a stress-rate based finite element algorithm for analysis 
of large deformations of inelastic bodies. It appears that the sole 
other analysis of large deformations of inelastic bodies by any 
similar algorithm is that presented by Atluri and Murakawa JH], in 


an invalid principle was presented by Levinson [ 1 6 J , and again 
by Zubov [17]. The failure of that principle is discussed by Dill 
[18]. 
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which necking of an elastic-plastic bar and postbifurcation analysis 
of a thin elastic-plastic plate was performed. The finite element 
algorithm used by those researchers was based on stress increments, 
rather than stress rates, and the motion of the elastic-plastic body 
was found by summation of increments. It was assumed that the 
accumulated error in this procedure could be kept small by minimizing 
'residual loads, 1 as may be done in elasticity. This procedure has a 
firm foundation for problems involving elastic bodies (whose deforma- 
tions were the subject of Murakawa's earlier research), but is of 
questionable validity when the body is not elastic. In their assess- 
ment of incremental solution methods for inelastic rate problems, 
Argyrls et al . [19] conclude that 

Inelastic rate processes are in general path-dependent; therefore, 
the drift (i.e. the accumulation of numerical integration errors) 
cannot be eliminated by residual load iteration, e.g. at the end 
of each time step. 

Moreover, when the body exhibits relaxation effects, this solution 
technique's numerical stability becomes extremely sensitive to the 
time step size. Hughes and Taylor [20] observe that the time steps 
required for stability in the explicit time stepping technique are 
much smaller than required for accuracy when only quasistatic deforma- 
tions are to be analyzed. 

A final objection to 'incremental' finite element formulations 
may be raised on the grounds that there always results an artificial 
coupling between the boundary value problem and the initial value 
problem. When dealing with 'flow law' type solids it is possible to 
treat the boundary value problem (for the rates) and the Initial value 
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problem (for the total stress and deformation) separately. Typically 
the boundary value problem for the rates is either precisely linear, 
or equivalent to a linear problem (without approximation). All of the 
nonlinearity falls into the initial value problem. Nonlinear initial 
value problems are perhaps the single type of nonlinear problem which 
we are best equipped to deal with numerically. In any case, we are 
better equipped to handle them than we are nonlinear boundary value 
problems. The incremental approach has the effect of actually trans- 
ferring that nonlinearity to the boundary value problem, where it is 
dealt with by residual load iterations. 

The objective of the present work is to develop a stress (-rate) 
based finite element algorithm for analysis of large quasistatic 
deformations of inelastic bodies. In doing so, we discard the notion 
of ‘increments' entirely. As a direct result, the boundary value 
problem and the initial value problem may be dealt with separately. 

The algorithm which results is applied to analyze large deformations 
of hypoelastic, hypoe 1 ast i c/pl as t i c, and hypoe 1 as t i c/vi scoplas t i c 
bodies . 

As is true in large deformation problems in general, the 
formulation of the boundary value problem and initial value problem 
is more complicated than in an infinitesimal deformation problem. 

The first part of this work is devoted to presenting, with reasonable 
completeness, the development of the p ro b lem. The finite element 
algorithm is not presented until Chapter VII. The initial value 
problem is presented in Chapter VIII. Example problems are presented 

in Chapter IX. We note from the outset that the finite element 
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algorithm is considerably more complicated, and involves more 
computation, th^n velocity based algorithms. However from the results 
it is clear that the improvement in accuracy over velocity based 
methods is substantial; so much so, that in view of the difficulties 
encountered in the application of velocity based methods to finite 
deformation problems, the present stress based approach appears to be 
the more efficient of the two. 
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CHAPTER I I 
KINEMATICS 

Kinematics is the study of deformations of bodies in space, 
without regard for the forces which cause them. Mathematically we 
represent natural space as a three dimensional Euclidean space E. 
Consider the motion of a body through space. The image of the body 
in E at the time t is the configuration C(t). As time passes, the 
configuration changes, and we say that the body deforms. In order to 
study aspects of deformations, we must be equipped to compare configu- 
rations assumed by the body at different times. We set out to equip 
ourselves thus. Let X be the position in E that was occupied by the 
material point X at the time T, and let x be the position of that same 
material point at the present time t. To indicate the dependence of 
x upon X, t, and t, we write 

x = X T (X,t) . (2.1) 

The function X T describes the deformation of the body relative to the 
configuration C(t) by 'tracking' each material point from its position 
in C(t). It is easily seen that at the time t * T 

X (X.t) I - x . ( 2 . 2 ) 

' 1 * t*T 


10 


ORIGINAL PAGE IS 
OF POOR QUALITY 


(2.3)* 


The gradient of X^ with respect to X 
F T (X , t ) = [V x X T (X,t)] T 

is called the deformation gradient, and is written for the 
determinant of F^: 

J T (X,t) = det F T (X,t) . (2.4) 

In view of (2.2), it is clear that 


F T (X,t) 


(2.5) 


t = T 


and 


J (X, t) I = +1 . 

1 1 t=T 


( 2 . 6 ) 


The time derivative of X T is the velocity function v^ : 




(2.7) 


The spatial velocity distribution is obtained by putting X^ 1 (x,t) 
for X in (2.7): 


'for a" 1 explanation of the special notations used in this work, 
see Appendix A. 
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( 2 . 8 ) 


-1 


v(x.O * y T (X x (* .t) »t) 


It Is clear that at the time t-T 


y„(X,t) I = v (x, t) . 
* t*T 


-T'- : 


We caution the reader by pointing out that v(x,t) and v (X,t) are 
entirely different functions. 

The velocity gradient L is defined by 


L(x, t) * [Vv(x,t) ] T , 


( 2 . 9 ) 


and it is clear that when t * T, 


F (X,t) I * L(x, t) . 
■" T " t“T 


The symmetric and skew-symmetric parts of L 


c - id + l t ) 


( 2 . 10 ) 


and 


w = i(L - L T ) 


( 2 . 11 ) 


have the physical significance of stretching and spin, and are thus 


12 


named. The trace of L has the physical significance of dilitation 


tr L * V • v . (2.12) 


The equat ion 


J T - J T (V* y) (2.13) 

is called Euler's expansion formula* in fluid mechanics. In view of 
(2.6), when t = T, 


J T (x,t) 


t=T 


V • v ( x,t) . 


We shall frequently write 3 for (7 • v) . 

Of course not any tensor field L is the gradient of a velocity 
field. The i ntegrab i 1 i ty condition (henceforth called compatibility 
equation) is 


Vx (l t > - 0 . 


(2.U) 


The ge-cfcl solution of the partial differential equation (2. 1*0 is 


L ' * 7y(x.t) . 


(2.15) 


Harris, lectures on fluid mechanics, Georgia Institute of 
Technology, Fall 1373. 
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Likewise, if 
tensor field 

V> 

then there i 
E 

an d 

OJ 


£ is a symmetric tensor field and w is a skew-symmet r i c 
, and (e - u) satisfies the compatibility equation 


< ( e - u) - 0 . 

-w **' * 


(2.16) 


s a twice differentiable vector field v for which 


*(Vv T + 7v) 


(2.17) 


i (Vv T - Vy) . 


(2.18) 
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CHAPTER I I I 


DYNAMICS 

Dynamics is "that branch of . . . mechanics which deals . . . 
with the action of forces on bodies in motion or at rest."* The 
fundamental laws of dynamics are called balance of linear mcneniur and 
balance of angular nonentum. They are expressed by the equations 

H - F = 0 (3-D 


and 


~ L - M = 0 (3.2) 

where H is the linear momentum of the body, L is the angular momentum, 
F is the aggregate force, and M is the aggregate moment. If no con- 
centrated or distributed couples are present, and if the reference 
frame is inertial, then H, L, F, and M are given in their classical 
forms by 

/ ft 

v(x,t)o(x,t)dV (3-3) 

V " 


American College Dictionary , Random House, 1970. 

V is the region in E occupied by the body at the present time. 
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(3. M 


where p Is 
intens I ty , 
body. We 
obtaining 


* 


L ■ f X X v(x,t)p(x,t)dV 
V ' 


F -/ b(x,t)p(x,t)dV + f TdS (3.5)* 

V 


M -/ xx b(x,t)p(x,t)dV +/ xxTdS (3-6) 

V ' " S ' 

the spatial mass density, b is the spatial body force 
and T is the true traction acting on the surface of the 
now introduce X T (X,t) for x in (3-3) through (3-6), thus 


H - / v (X,t)p (X)dV 

V T " 

L ■ / X T (? x y T (X,t)p (X)dV 

■ I# “■ * L 


(3-7)** 


(3.8) 


F - / b(X,t)p(X)dV +/ TdS (3.9)*** 

V t- $t 

M =/ X_(X,t) xb (X,t)p (X)dV +/ X T (X,t)xJ dS. (3-10) 

• , T i L T • T 



S is the surface of V. . 

is the region in £ occupied by the body at time T. 

S is the surface of V . 

T T 
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The terms P T (X ) , b T (X,t), and J are defined by 


P x (X) = p(X,T) 


b x (X,t) * b(X T (X,t) ,t) 

and 


J T = J(dS t /dS x ) . 

We shall consistently refer to b^ as the nominal body force 
and to as the nominal traction. Inserting (3-7) through (3-10) 
into the dynamical equations (3-1) and (3-2), one may deduce by the 
usual arguments that 


T + T = 0 

-I -T 


(3-11)* 


It = 0 T - t T 


(3-12)** 


7 • t_ + o_b. 


= p v 

T-T 


(3-13) 


r _ ♦ t _ - (F_ • t_) ■ 0 

~ T -T ■**- T * s ' v *■* 


(3-iM 


the + and - indicates that the two tractions act on opposing 
sides of a surface with a well defined normal. 

"V(X) - n(#,t) 
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called traction reciprocity, the stress principle, and the local forms 

of linear momentum balance and angular momentum balance, respectively. 

The stress t Is called the nominal stress.* If the reference con- 
-T 

figuration and present configurations coincide (t*x), then It is 
eas i ly seen that 



I(X.T) 


and 


t (x,t) j - T (X, t) . 

t=T 

The familiar dynamic equations are thus recovered: 

J + + j" = 0 


J » n • T 

V • T + pb * pv 


(3.15) 

(3-16) 

(3.17) 

( 3 . 18 ) 


called the Lagrange stress by Prager [21], the First Piola- 

Ki rchhoff stress by Truesdel 1 and Noll [22J; other def initions result 

if T ** t *n instead of (3.12). 

-T ~T -T 
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where T Is the true stress.* 

According to their definitions, the nominal tractions and 

obey 


V S T ■ - T ? dS ? 


(3.19) 


Using the stress principle (3-12) and the kinematical relation 


n ; dS^ = J X ^)0 T * <5)dS x 


( 3 - 20 ) 


we obtain from ( 3 - 19 ) the relation between the nominal stresses t_ 
and t^: 


t_(t) * j t (Of^(;) • t ; (t) 


(3.21) 


When one of the configurations C(t) or C(c) coincides with the present 
configuration, say £ -*• t , then we obtain 


t T (t) ■ J x (t)F’ 1 (t) • T ( t ) , 


( 3 - 22 ) 


relati^-7 fh? t n.i B stress T to the nominal stress t^ . Of course, 
equations ( 3 - 19 ) through (3.22) are for application at a material 
point (as opposed to a spatial point). 


»V 

also called the Cauchy stress. 
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We now set ( 3 - 1 0 through (3.1*0 for traction rates and stress 
rates.* By stress rate we mean the actual time derivative of stress. 
This distinction is necessary because in some instances (not in this 
work) the word 'rate* is used interchangeably with the word 
'increment. 1 In this work a superposed dot indicates a material 
derivative; that Is, the time rate at a material point (as opposed to 
a spatial point). 

The rate forms of (3-11) through (3-1*0 are found by ordinary 
differentiation. We thus obtain 


• + * - 

Jt + - T X - -° 


T * n • t 
-T -T ~T 


V X ' Jt + ' °T- V T 


( F • t +F*t)-(F*t_ + F*t) =0 
~T ~ T ~T ~T ~T -T ~T ~T 


From equations (3-21) and (3-22) we obtain 


■ - 1 * 

It * J t!t * ! 


(3.23) 

(3-2*0 

(3-25) 

( 3 . 26 ) 


(3.27) 


this discussion is unrelated to the controversy surrounding 
the definition of an 'appropriate stress rate' for constitutive 
theory . 
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and 


t - - (e + oj) • T + Jt + t . 


For the stress rate t, equations (3-23) through (3-26) become 


• 4 > • — 

It + It ■ 2 


T = n • t 


• • 


V • t + pb = pv 


(E+U))*T + t- t - T*(e-Oj) =0 


( 3 . 28 ) 


(3-29) 
(3-30) 
(3.30 
(3- 32 ) 


Equation (3-27) will prove essential in the eventual integration of 
the stress rate. Equation (3-28) is important in the reformulation 
of the constitutive equations. It is clear that at the time t = x 


t (X,t) I - t(X,T) . 
t — T 


We now discuss the difference between the nominal traction 
rate and the true traction rate T. Differentiation of (3-18) yields 


T *n*i + n*T 


(3.33) 


The rate of change of the normal vector n is 
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(3.3*) 


a 

n = n • [- (e + oi) + (n • e • n) I ] , 

Elimination of n from (3.33) yields 

T = n • (-(c + u) • T + (n • e • n)x + xj , (3 - 35 ) 

From (3-28) [- (g + u) * t + tJ = (t - Jx) , and since n • t « T , and 
n*x = J, we finally arrive at 

I * J t - (j - n • e • n ) T , (3-36) 

Though (3-36) was arrived at by manipulation, its physical significance 
is easily found out. Consider a force P acting uniformly over a flat 
surface S of area A. Then the nominal traction rate T is given by 

T t - (P/A) , 

and the true traction rate T by 

T - (P/A) - (A/A) (P/A) . 

However, since (P/A) ■ T and (P/A) m J, this last equation may be 
written as 

T - T - (A/A) J . 
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Glancing back at (3-36), it is now clear that the difference between 
the true traction rate and the nominal traction rate is due to strain- 
ing of the surface. In applications a hydrostatic pressure may be 
prescribed on a straining surface. In such a case the formula (3-36) 
must be used to relate the nominal traction rate to the true traction 
rate . 

Some other practical information can be drawn from equation 
(3-36). Suppose part of the surface of a body is to be t ract ion- free 
during a deformation; that is, T(t)*0. Then T * 0 is the appropriate 
boundary condition for the nominal traction rate. The boundary con- 
dition T * 0 and equation (3-36) yield 

T + (j - n • e • n)T * 0 ; T(0) - 0 . . (3-37) 

“ * -V “ ” “ ” ” 

This is an i n i t i al value prob lem for T(t). The solution is T( t ) « 0 . 
However, in a numerical integration of the boundary value problem we 
should expect errors in J to be amplified or attenuated as the surface 
contracts or expands* (that is, as the coefficient of J above is 
negative or positive). This example is simple, yet suggestive of the 
types of problems associated with traction boundary conditions on 
surfaces which experience large strains. 

For norational convenience we define the Ki rchhoff stress 0 T as 


the stability of solutions of the traction reciprocity 
equation (3.29) depends on the surface stretching in the same 
manner. 
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( 3 - 38 ) 


?T “ J Tl * 


It Is clear that at the time t-T 


a T (X,t 


t*T 


t(X,t) . 


( 3 . 39 ) 


The rate .of the Kirchhoff stress Is related to the rate of the true 
stress by 


a T « J t t + J t t 


( 3 -^ 0 ) 


and it is clear that when the reference configuration and present 
configuration coincide 


• • • 

a - Jt + T 


( 3 - 41 ) 


From the last two equations and Euler's expansion formula we find the 

• • 

relation between and a to be 


o T * J T a , 


( 3 . 42 ) 


and it is clear that when t»T 


't* 8 ! 


a{ t ) . 


( 3 . 43 ) 
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The relationship between the Kirchhoff stress rate and the nominal 
stress rate may be drawn from (3-28) and (3-^1) as 

t * -(e + w) • t + o . (3.^) 

This equation shall be used in the reformulation of the constitutive 
equat ions . 

The final topic in this chapter is construction of the general 
solution of the linear momentum balance (3-25). That equation is of 
the form 


D l f l + °2 f 2 * 0 (3 ‘ 45) 

where D 1 = V y * , D 2 » f 1 * t T> and f 2 = P t " -t^ ’ Just 35 

for differentials in'the plane, the general solution of (3-^5) is 
( forma 1 1 y ) 


fj ■ D 2 <}> + k j ; f 2 * + k 2 ; (3*^6) 

where 4> is arbitrary and k. is any function which satisfies 
D.k. = 0 (no 5 urn) . 

i i 

It follows immediately that the general solution of the rate form of 
the linear momentum balance is 
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(3-^7) 


t T (X,t) - $ T (X,t) + V x x # (X,t) 

P T (X)(b T (X,t) - y T (X, t) ) - -7 X • ^ T (X,t) + * T (X) ( 3 .A 8 ) 

where 3^, and y^ are arbitrary functions of the indicated arguments. 

From our experience with natural bodies, we often expect a 
gradually applied statically balanced load to produce a 'gradual' 
deformation of a solid body. The obvious counterexamples are snap- 
through phenomena, but it is the very failure of the gradually applied 
load to produce a gradual deformation which leads us to call such 
behavior 'unstable.' The general mathematical model for the motion 
and deformation of bodies is not so well understood that one can 
ascertain the stability of the motion of a body (subjected to gradually 
applied loads) without actually determining that motion and 'inspecting' 
it. Of course this is not practical, so it has become common to 
anticipate stability (for gradually applied loads). 

if we assume that, for sufficiently gradual application of 
loads, the ensuing motion of, a body is stable, then it is reasonable 
to ignore the inertial terms H and X in the dynamical principles (3-1) 
and (3-2). This is known as the quasistatic assumption (or less 
accurately, the quasistatic approximation). In such a case the 
general solution of the linear momentum balance may be simplified.* 

We write t as the sum of a homogeneous and a particular solution 


* 

only 


the quasistatic solution of ( 3 * 31 ) 


is recorded here. 
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(3.*9) 


• *o *b 
t - t + t 


•o . 


where t is the solution of the homogeneous equation V*t*0, 


t°(x, t) =* Vx *(x, t) 


and t 


is any particular solution of 


(3-50) 


V* t b (x,t) = -p(x,t)|~ b (x, t ) + v(x,t) »7b(x,t)j . (3-51) 


• b 

A particular solution t may be constructed in cartesian coordinates 
using indefinite integrals 


t 


b 

' j 



+ v 



(no sum) . 


(3-52) 


• ^ 

Notice that t will depend upon v unless v*Vb = 0. If the only body 

force is due to gravity then one sets Vb = 0. However, if D'Alembert's 

* b * 

principle is used, then t will depend upon v»Vv.* 

The angular momentum balance for t (3-32) involves the stretch- 
ing and spin, as well as t. For this reason we cannot form the general 
solution of that equation by judicious choice of stress functions, as 
may be done in the case of the true stress. The function $ is there- 
fore called a first order stress function. 


in a steadily spinning disk v • Vv ■ 0 , since the radial velocity 
vanishes and the acceleration is entirely radial. 
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CHAPTER IV 


CONSTITUTIVE EQUATIONS 
Int roduct ion 

A material is characterized by its response to loads and 
deformations. The mathematical model for the load- response behavior 
of a material is called its constitutive equation. A constitutive 
equation of the form 

L * B ( t , T ) 

would enable us to set a boundary value problem involving the stress 
rate t alone.* The possibility of finding such a constitutive 
equation is now discussed. Consider a deformation for which L, t, and 
T are the observed velocity gradient, nominal stress rate, and true 
stress. If a rigid motion 

x' ■ Q(t) • (x - x ) + x, (t) 

is superposed upon this deformation, then the apparent velocity 
gradient, nominal stress rate, and true stress are L', t' , and T 1 


( 4.21 


(4.1) 


this form has appeared in the literature; see Hill [23]* 
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(distinguished by a prime). The tensors L', t 1 , and x' are related to 
the tensors L, t, and T by the formulas 

L 1 * • L • Q T + Q • Q T (4.3) 

•W *S» <V 

(4.4) 

(4.5) 


• . • T • T 

t ' * Q • t • Q. + Q • x • Q 

x ' = Q • x • Q T . 


The principle of material frame indifference states that the mechanical 
behavior of the material is indifferent to rigid motions such as 
(4.2). Therefore, the constitutive function B must satisfy 


L' * B(t ' ,x 1 ) 


(4.6) 


for the same function B. To establish whether or not B satisfies 
frame indifference, we use (4.3) through (4.5) to eliminate the primed 
quantities from (4.6), thus obtaining 

Q • B ( t , x ) *Q T + Q • Q T = B ( Q • t * Q T + Q • x • Q T , Q * x • Q T ) , 

(4.7) 

which must be satisfied by all orthogonal Q and all skew symmetric 
Q*Q at everv moment of time. In particular, it must be satisfied 
when 
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i ; 


• *T 
Q « -Q . 


(4.8) 


Then (4.7) serves to restrict the form of B: 


B(t ,t) * B(t + T • Q T ,t) - Q 


(4.9) 


for arbitrary skew symmetric Q. The right hand side of (4.9) depends 

upon Q but the left does not. Thus (4.9) can be satisfied (for 

• • 

arbitrary Q) only if Q cancels out on the right hand side of (4.9). 
This is not possible in general; that is, no functional form for B 

assures that Q. cancels.* 

«»■ 

This difficulty may be overcome by postulating a new constitu- 

.. • 

tive function B which depends not only upon t and T, but also upon 

u. ** Then we easily replace (4.9) by 


B(t ,t,'cj) * B(t + T • Q ,T,us+ Q) - Q , 


(4.10) 


to be satisfied for arbitrary skew symmetric Q. Again, this is only 
possible if Q. cancels out on the right hand side. Therefore B must 
be of the form 


B(a,b,c) ■ F(a+b • c,b) + c ; 


(4.11) 


A 

further light is shed on this problem by Ogden [24], 

* or for^that matter, any tensor T which transforms by the rule 
T' - Q • T • qJ + Q • Q t . ~ 
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i n other words , 


L = F(t + T*io,l)+u). (A. 12) 

*** %* ^ *** 

Defining a stress rate p as 

p * t +T • to (*4.13) 

and recalling that L-lo=c, we replace ( 14 . 12 ) by 

e = F(p,x) . {b. ]k) 

• ^ 

The stress rate p is called the Lure stress rate in this worl; * Frame 
indifference requires that F satisfy, for all orthogonal Q, 

q- ? ( e* 3 } * * p • S T ’S’ j • s T ) • (*»• 15) 

A tensor function possessing this property is called isotropic. From 
the theory of matrix polynomials [25], we know that if F is a sym- 
metric isotropic function, at most affine** with respect to p, then F 


* • 

p has the physical significance of a 'corotational nominal 
stress; Tt has been used by Wunderlich and Obrecht [9] for beams. 

an affine function is composed of linear and inhomogeneous 
parts; for example, if f(x)»ax + b, then f is affine w.r.t. x so long 
as a and b are independent of x. 
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may be represented as a linear combination of the six tensors 
T 1 s 

«w 

Hr • T' + T' • r) i(r • s + s • r) , 

• • • T \ 

where r= i(p + p ) and s*T' • t 1 . However a more useful representation 
of F can be found which depends upon the traditional form of the 
constitutive equation. Thus, we defer further discussion of F until 
after the traditional form has been presented. 

Traditionally one assumes a constitutive equation of the form 

a = £(e, u),t) . (^ . 16 ) 

It is a simple application of the frame indifference principle to show 
that E must be of the form 

(*.17) 


(*.18) 


where C is a symmetric isotropic function. For convenience we define 
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£(a,b,c) * C(a,c) + b • c - c • b 


so we replace (**.16) by 


c-u'T + T'U )* 1 C(e,t) , 


I 


r 


IIUAI H IdlllliiU i I i ,|| II ill III IN. i || 



the corotational stress rate a* by 


a* 


a - uj • t + t • to 


( 4 . 19 )* 


If C is at most affine with respect to £, then the general form for 
(4.18) is 


a* * V: £ + Z 


(4.20) 


where the fourth order tensor V and the second order symmetric tensor 

Z depend upon the stress. The components of V and Z, in cartesian 
— * — 

coordinates may be set in the forms 


V . . - X 1 1 (6. . 6. , ) + X ,z (6..t' ) + X I3 (6..s. .) (4.21) 

ijkl ijkl ijkl ijkl 


12 


1 3 / 


+ x 2 ' (T ij 6 ki> +x 22 h i'j\i> 


+ X^(s..6 .) + X^ 2 (s..t‘ ) + X^(s..s .) 

ijkl ijkl ijkl 


♦ Jl;,({ ikV + ' , 2 (s ik T ij + 1 ik { ij ) * lj 3 ( 5 ik s ij ts ik s ij ) 


the essential property of this stress rate, or any other stress 
rate described as 'corotational 1 in this work, is that it is absolutely 
insensitive to rigid spin of the body. In frames in which the spin 
vanishes, such stress rates always reduce to ordinary time derivatives 
of stress. 
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and 


E ij “ ^ij + r,2T j j + n Js |j • (4.22) 

In (4.21) and (4.22) T . is the stress devi a tor, 4. r * T ' T.' . , and 

i J i j ik kj 

i I i i _ _ . . 

A , y , and n are functions of the invariants of t. If A IJ «A-* f 

then •' V k j .j, and we describe V as symmetric. Equations, (4 . 20) 

through (A. 22) include, as special cases, the theories of h'ypoelastic 

and hypdelastic/vlscoplast ic bodies , and with sTight modification, 

plastic bodies. These shall be discussed individually at the end of 

this chapter. 

We now construct from (4.20) a representat ion for F (4.14). 
Recalling from ( 3 - ^« ^ ) that the nominal stress rate and the Kirchhoff 
stress rate are related aS 

t - [c - (e + w) • t] , 


it is easily seen that the relationship between the Lure stress rate 
and the Kirchhoff stress rate is 

p - t + t*w - [a - (e + y) • f J + t*w • 

this yields, using the definition of the cdrotit ionaf stress rate, 

p * g* - e • t . (4. 2 3) 
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The symmetric and skew symmetric parts of p are given by 


Hp + p T ) - a* - He • t + t • e) 


(4.24) 


and 


Hp - p T ) 


■j(e 


e) . 


.25) 


According to angular momentum balance, the latter of these (4.25) must 
be an identity. Using (A. 20) to eliminate a* from (A. 24) yields 


= w :e + Z 


(4.26) 


where W is defined by 


W 


i j k 1 ^ i j k 1 


J(i., T.' . +l!. 5, .) 

i k 1 j i k 1 j 


T T nr>m (<S i k 6 ! j ’ ’ f' 1 - 2 ?) 


Note that if V is symmetric then W is symmetric also. 
= * 

The stress rate r 


r = i(D + p T ) * •'(t + T*u-U) , T + t" r ) 


(4.28) 


is called the symmetrized Lure stress rate in this work. It is the 
time rat® of a symmetric nominal stress used recently by de Veubeke 
[26], and called by him the Jaumann stress. Since the corotational 
stress rats f often called the Zaremba- Jaumann rigid-body stress 
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rate, different from the rate of the stress used by de Veubeke (r), it 
would be likely to cause confusion if the name Jaumann were assigned to 
any stress rate In this work.* 

If W is invertible, then the constitutive equation (*.26) may be 

ae 

wr I tten 


£ =* W _1 ; (r - I) (*.29) 

and we see that the function F (*.l*) corresponding to (*.20) is 
represented by 

F(p,t) - W 1 : [3 (p + p T ) 

where W is defined by (*.27). 

In applications one is usually given the constitutive equation 
in the traditional form (A. 20), (*.21), and (*.22). To obtain the 
form (*.30), one must construct W according to (*.27) and then invert 
W, if possible. This is a major undertaking from a computational point 

as 

of view, for W will generally be different at each point of a stressed 
body, even when V is constant. Therefore special attention is given to 

A- 

practical methods for construction of W For plane problems W ' can 

* % 

be found in closed form. For general problems in which V 1 is known 

ft 

a truncated power series for V * is often of acceptable accuracy. The 

ft - ■ 

details of these two special cases are discussed in Appendix B. 
see At luri [153- 


- n. 


(*•30) 
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Hypoelastic Bodies 


If E*0 in equation (4. 20), and the coefficients X'-^ and y' are 
analytic functions of the stress Invariants, at least in some open 
region of stress space (containing the origin), then the body is said 
to behave hypoelast ical ly in that region [26]. No relaxation phenomena 
will be exhibited by the body so long as E»0. The state of stress at 
a material point will depend upon the deformation history, but not upon 
the speed of the deformation. This makes it possible for us to use any 
monoton ica 1 ly increasing parameter, such as a characteristic displace- 
ment, in place of natural time when studying the deformations of 
hypoelastic bodies. 

The first approximation to the behavior of any hypoelastic body 
for small deformations from the stress-free state is obtained by setting 
T = 0 in ( A . 2 1 ) . Then V may be written as 


v ijki ■ x(s ijV * 2u(5 ,kV ■ 


(Ml) 


where X (Lame's constant) and y (the shear modulus) are defined 


X = X 


11 


y * y 


T=0 


T=0 


Referring to (4.27), we construct W as 


W... . - X(6. .6. .) + (2y-i t )(«.. 6..) 
ijkl ijkl 3 mm iklj 


(4.32) 


' i(s ik T ij + T ikV ■ 
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The constitutive equation given by (4.31) must be regarded as 
the simplest possible model of solid behavior. It is doubtful that any 
real material behaves as predicted by (4.31) when deformations are 
finite. Its popularity stems not only from Its simplicity, but from 
the success it has enjoyed in infinitesimal strain theories, and the 
lack of physical data needed to use more sophisticated models. Never- 
theless, there is a point beyond which any semblance of real material 
behavior vanishes from (4.31). Consider rectilinear shearing from a 
stress-free state: 

1 ? ? 3 

v * kx ; v ■ v » 0 ; T “0. 

~ t*0 

The constitutive equation (4.31) predicts the following stresses [27]: 

T ]2 * p s in (kt) ; Tjj » -T 22 » V ( 1 ” COS (kt) ) . 

Until the time that kt ■ i^r, the shear stress t 12 increases. After 
that time Tj 2 decreases, even though the shear strain continues to 
increase. In this particular problem we may take (kt) to be a measure 
of the strain; beyond the strain kt ■ (4.31) fails to provide an 

acceptable model of any material's behavior. 

The rectilinear shearing example above is such a simple 
deformation that we were able to rely on intuition alone In deciding 
the limit of applicability of the equation (4.31). For more compli- 
cated deformations (and/or more complicated materials) our intuition 
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is not so strong. The notion of 'realistic material behavior' needs 
to be set down in a form which (minimally) allows us to distinguish 
between 'realistic' and 'unrealistic' deformation processes. This is 
the role played by thermodynamics in modern continuum mechanics [28]. 

Yield Surfaces 

It is well known that metals exhibit behavior which is more or 
less elastic so long as strains are small. In other words, when 
moderate loads are imposed, then removed, the elastic body returns to 
its original shape. If more severe loads are applied to the body, then 
removed, some permanent distortion of the material may occur. The 
mechanisms of inelastic behavior may become active very suddenly, as in 
mild steel, or only gradually, as for lead or copper [29]. In any 
case, one way of idealizing the transition is to introduce a yield 
surface to stress space, inside of which the mechanism of inelasticity 
is dormant, on or outside of which the mechanism is active. The yield 
surface may change during the deformation as dictated by the change in 
the yield behavior of the real material being modeled. After the 
initial yield of the material, the surface is usually called the load- 
ing surface. 

The von Mises' yield criterion is the most common of those found 
in the engineering literature. The equation of von Mises' surface is 
expressed by either of 
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J 0 VJT r - o’ 1 (k) - Oi 

J 0 '^~i 7 ’ °I 2(k) ■ 0 ■ 

(1 1 ■ t 1 : t 1 ) 

^ ■»» 

1112 

where a and a are material functions called the effective uniaxial 
e e 

stress and effective shear stress, respectively. The parameter k is a 

monotonical ly increasing invariant of the local plastic deformation 

history. The function c^(k) must be determined from uniaxial test 

data, and the function a g (k) must be determined from plane stress or 

strain rectilinear shearing test data. The scalar J q is the volumetric 

strain relative to the stress-free state. It must be included in the 

first of (1*. 33) because in uniaxial tests, the total load P and length 

1 are measurable. The initial dimensions of the specimen are known, 

so the measurable stress is : 

o 

a 11 .jt"> (A1/A 1 ) (P/A) - P(l/A 1 ) . (4.3*0 

o o o o o o 

Similar arguments may be given for inclusion 6f J q in the second of 
(4.33)- The usual Mises 1 criterion is recovered by dividing (4.33) 
through by J q . In practice would usually be expressed as a 
function of the mean stress: 

J = (3a + 2u)/(3X + 2y - t: I) . (4.35) 
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Of 


course a"and a' 2 
e e 


cannot be 


independent ; 


in view of (4.33) 


they must be related as 


o' 1 (k) - O l2 (k) . 
e e 


(4.36) 


If physical dat£ cannot be reconciled to this condition, then the 
body in question does not admit von Hises' representation of the load- 
ing surface. 

One might suspect that von Hises 1 criterion Is no more than an 

extreme idealization. However, according to Hill:* 

experimental data for mild steel . . . suggests that the yield 
locus changes over from a hexagon to a circle with progressive 
cold work. However, for other steels, and for copper and aluminum, 
von Hises' criterion appears to fit the data equally wel 1 no mat- 
ter what degree of prestrain. 

So, while it is indeed an , i deal i zat ion , its agreement with physical 
data for some metals is quite good. Nevertheless, the small dif- 
ferences may be important. According to Chr i stof fersen and Hutchinson 
[ 31 ], who have proposed a class of 'corner theories' of plasticity, 

it Is generally agreed that the simplest flow theory built upon 
the assumption of isotropic hardening using the Hises yield sur- 
face underestimates certain crucial plastic strain components in 
a non-proportional stress history such as encountered in buckling 
or necking. 

In contrast to a smooth yield surface characterization, a corner 
theuiy will niosi likely overestimate certain components of the 
plastic strain increments in the vicinity of an abrupt change from 
proportional loading. 


*H ill [30], p. 24. 
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As a consequence, "calculations based on a flow theory with a smooth 
yield surface give necking-type bifurcations at strain levels which 
far exceed realistic values" [ 31 ], while calculations based on corner 
theory can be expected to provide conservative estimates of the 
bifurcation strain. Finally, according to Tvergaard et al . 132], who 
have presented a numerical study of flow localization In t£e plane 
tensi le test, 

Analyses carried out within a theoretical framework . . . reveal 
that the classical elastic-plastic solid with a smooth yield sur- 
face is quite resistant to the localization of deformation Into a 
shear band. 

The resistance of the classical elastic-plastic body to intense local 
deformation may even be sensitive to the yield criterion's independence 
from the mean stress (such as that induced by a distribution of voids). 
Thus, the circumstances under which von Mises' criterion can be 
expected to provide reasonable results have not been resolved. 

We observe that, as a rule, smooth deformations are found when 
a smooth yield criterion is used. This may be an over-simplification, 
but It is borne out by the weight of solved problems in the literature. 
Regarding finite elements, and the goals of the present work, this is 
important because we may be reasonably confident that intense local 
deformations, which would require a much finer finite element mesh 
than has been used, do not occur. In view of the discussions of 
Tvergaard et al. [32], it would seem less safe to make such an 
assumption if the most common alternative to von Mises* criterion, the 
Tresca 'maximum shear stress criterion,' were used, since that surface 
has vertices. In this work, von Mises' criterion is used exclusively. 
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The parameter k Is usually chosen as either the 'plastic work' 

W p or some variant of Odqulst's 'plastic strain' q, defined by 
Kachanov [29] as 

w p « f *t : e p dt ; q - / * Jat p : e p dt . (4. 37) 

Notice that both of these Increase monotonical ly during plastic 
deformation. In uniaxial extension 


W P 





or in rectilinear shear (plane strain) 


W P 


„ 12 P 
2T ^ 



In any case we can find a functional dependence between W p and q 

3 

since 


dW P /dq ■ (t : E P )/(aE P : £ P )* > 0 , 

3 


so the choice of parameter is largely a matter of convenience. 

Hypoelastic/Plastic Bodies 

We give a brief sketch of the theory of hypoe last i c/plast i c 
bodies. The reader is referred to Hill [30], Prager [33]» and 
Kachanov [23] for extended treatments of the theory. 
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As was discussed In the main text of this chapter, the 
mechanical behavior of an Isotropic body Is completely known once 
material functions X^ t y', and r|‘ are known (see equations . 20 
through A. 22). However, no series of physical tests can be devised that 
would determine these material functions completely. In the face of 
such indeterminancy, it is customary to choose A J , y , and n’ which 
produce the simplest constitutive equation capable of explaining the 
behavior which has been observed.* The classical theory of plasticity 
for polycrystalline metals is just such a theory. 

Suppose the plastic stretching e p depends upon the stress and 
stress rate as 

e p « v' : a* + Z' 

where V 1 and Z' are of the forms (4.21) and (4.22). Then from the 
apparent absence of relaxation phenomena, at least for quasTstatlc 
cold-working of the material, one concludes that Z 1 must vanish. The 
physically observed incompressibility of plastic deformation places 
the following restriction on V ' : 

I : V' - 0 . (**.38) 


It is entirely possible that (**.20)- (**.22) and the physical 
data cannot be reconciled; in such a case one must discard (4,20)- 
(*♦. 22 ). 
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From a series of uniaxial tests in which the Ki rchhof f stress o'* is 
found as a function of the logarithmic strain 

- '"('/'o’ 

* 1 1 

we establish the following relation between o and 

o 

o'Ve,, = {J E ( k ) } when a" = a* * (k) 
oil o t o e 

and o* * E, , >0 

o 1 1 

c'Ve.. ■ {J E ( k) } otherwise 
oil o 

which is used to further restrict V 1 . The material 
general representation for V 1 cannot be determined uniquely from the 
restrictions (4.38) and (4. 39). However, a very simple solution may 
be found as follows: we assume that 

V 1 - A 22 t't' . (4.40) 

mt — 

22 

This form satisfies (b. 38 ). Now we define A so that (b. 39) Is 
satisfied. The plastic stretching is given by 

A 22 t' (t' : a*) » e p . (4.41) 

The scalar product t' : e* 3 is therefore 


1 1 


(4.39) 


functior.5 in the 
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A 2 2 ( T ' : T ' ) (t 1 : o*) - T* : e P . 

** mr 

Putting I ' - (t ‘ c T ' ) , I ' *■ J 2 I ' , and W P ■ T : £ P , from this last 
equation we find that 


A 22 * (2/ 1 * ) (dW p /d I * ) . (4.1*2) 

The function (dW P /dlj.) may be determined from the uniaxial test data 
(4.39) as follows: 


dvr 


11 p 

T £, 


1 1 


dt 


1 4 11*11 

d f t * y t 'o"dt ; 


so 


(dW p / d l ^ ) - | (e p ,/o n ) . (4.43) 

Recalling that O o *J Q a, we get from the test data 

(e p /o 11 ) - J Ke-e^/a’ 1 ] - [(E.(k)}" 1 - (E(k) }” 1 3 

I J O O t 

= (1/h) . (4.44) 

On replacing o^ 1 by J q Jj l~ , and recalling (4.42) and (4.43), we 
2/ 

finally obtain A as 
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A 22 - (j l') 1 (^-)[{E t (k) > 1 - { E ( k ) } '] for loading 

(i f JqVP ■ (k) and t : e ^ oj 

* 0 otherwise ( 4 . 45 ) 

or, in terms of the function h, defined in ( 4 . 44 ) above, 

A 22 ■ (y *'^) ^ ( j~) for loading ( 4 . 46 ) 

■ 0 otherwise. 


Values of the function E (•) are called the tangent modulus. 

n o 

Values of E(*) are the instantaneous Young's modulus. The function A 4 

is multivalued on the loading surface; this is the sole source of non- 

linearity in plasticity theory. If the dependence of A upon k were 
22 

removed then A would be single valued on the yield surface. The 
resulting hypothetical material is encountered in the context of 
uniqueness theory, and is called the 'linear comparison solid' there. 

The complete constitutive equation is obtained when we consider 
the total stretching* to be the sum of elastic and plastic parts: 




V 1 : o* 




< fi lk«lj> ' I 


ij-kP + A 


2 v T * 

T 1 j T kl 


A 

?.e., Lhe svnm'cLiic part of the velocity gradient. 


( 4 . * 7 ) 
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i 


Notice that V ^ thus defined is symmetric; that is 

St 

v ” 1 ■ v“ i 

Ijkl v klij ’ 

22 

Since A is multivalued on the yield surface (4.47) is not linear^ 

We note the agreement between the present equations and those of 
McMeeking and Rice [34]. 

Hypoe lastic/V is co plastic Bodies 
The theory of vi scopl ast i c E ty is similar to classical 
plasticity theory in the sense that it provides the simplest constitu- 
tive equation capable of explaining the uniaxial test data of certain 
materials. The difference is that the materials being modeled exhibit 
creep and relaxation. The viscoplastic model common in the literature 
Is compatible with 

isotropic yield behavior 

plastic incompressibility 

uniaxial creep and relaxation tests. 

The representation of the viscoplastic stretching c Vp found In the 
monograph of Perzyna 1 35 3 is of the form 


E 


vp 



(4.48) 


where y is a function of the stress invariants and perhaps the 

deformation history. It may be determined from uniaxial test data In 

22 

basically the same manner as A was determined for a plastic body. 
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Suppose that from uniaxial tests we establish 




11 


when a 11 > a' 1 (k) 
o e 


- 0 


otherwise. 


(4. k3) 


For uniaxial tests, the general equation (4.48) reduces to 


eiT ■ \ yc4 Vj when a* 1 > (k) . 

1 1 3 o o o e 


Using the result of the uniaxial tests, we find y as 

e vp 

Y * \ J -44* - \ J i»(cr ] ' ,k)/a' 1 when o' 1 >a ^ (k) 

c o _l 1 z o o o o e 

a 

o 


* 0 otherwise . 


( 4 . 50 ) 


(4.51) 


The complete constitutive equation is obtained when we consider the 
total stretching to be the sum of elastic and viscoplastic parts: 


- e e ♦ E v|> - V ’ 1 : i* * 


YT 


where 


V i jkl “ ( € ) (<S ik 6 lj* " E (<S i j 6 kl ^ 


( 4 . 52 ) 


(4.53) 
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and y is defined in (it. 51). Notice that V ' is symmetric. Only this 

m 

simple model of viscoplasticity is used in the examples accompanying 
this work. 


50 


CHAPTER V 


BOUNDARY VALUE PROBLEMS 

In the preceding sections we have treated kinematics, dynamics, 
and material behavior as distinct subjects. We obtained equations of 
compatibility (2.16), linear momentum balance (3*30, angular momentum 
balance ( 3 * 32 ), and the constitutive equations (4.20) and (4.26). 
Presently we regard these equations as a system of coupled partial 
differential equations. It is worthy of special note that all of 
these equations are linear, with the exception of the constitutive 
equation in the case of a plastic body. For ease of reference we col- 
lect these equations below: 


V x (e - to) * 0 in V; s • (-£ + u> + 7v) * 0 on S ; (5 - 1 ) 

V* t + pb « 0 in V ; n • t ** J t °n S ; (5.2) 

iI(e + <i})‘T + t-t T -T«(e-u)]*0 In V; (5-3) 

t-Vre-e^T-T'w + E; (5.4) 

^ 

c * W 1 : [i(t + T • u) - ui • r + t T ) - E] ; (5.5) 

y ■ v on S v ; (5.6) 
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(5-7) 



These equations have been discussed Individually in the preceding 
chapters, with the exception of the second of (5.1), and (5. A). The 
constitutive equation (5.M follows directly from (*4.20) and 

t ■ a* - e * x - t • w . 

*» *v ^ ^ 

In the second of (5-1) s is an arbitrary surface tangent; the equation 
is essentially a counterpart to the stress principle, as it relates the 
(two-dimensional) surface-velocity gradient to the (three-dimensional) 
velocity gradient. We shall Imply only symmetric tensors by writing 

e: e - = 0 (5-8) 

V "m 

and oniy skew symmetric tensors by writing 

a): u + w T » 0 . (5.9) 

This is rather trivial, but nonetheless necessary to complete the 
system of equations (5.1) through (5-7). For the sake of clarity we 
shall often refer to (5-1) as 'compatibility,' (5.2) and 'LMB' (linear 
momentum Da lance ) , (5-3) as 'AMB' (angular momentum balance) , (5.6) as 
'VBC' (velocity boundary condition), and (5-7) as 'TBC' (traction 
boundary condition). We call the system of equations (5.1) through 
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(5.9) 'the general boundary value problem.' 

As discussed, the general solutions of compatibility (5.1) and 
LMB (5.2) are given by 

( 5 . 10 ) 


. 

pb / 

respectively. The remaining equations of the general boundary value 
problem (5-3) through (5.9) are algebraic in character. Derived 
boundary value problems are those obtained from the general boundary 
value problem through use of (5.10) or (5.11). Usually one uses (5.10) 

to eliminate e and u> as variables from the constitutive equation (5.**). 

*%» 

so that £, cl), and t are all determined by v. Owing to the special 
structure of the constitutive equation (5.**), AMB Is satisfied 
implicitly for all such e, w, and t. Upon elimination of t from LMB, 
one obtains a single second order partial differential equation In v. 

We call this the First boundary value problem. Alternatively, one may 

use (5*11) to eliminate t from the constitutive equation (5.5), so that 

• . 

£ is determined by <f> and u). Subsequent introduction of 4 (via t and e) 


£ - }(Vv' + Vv) 


oj . $ (Vv - Vv) 


and 


• *o *b 
t - t + t 


t“VX$; V»t 


53 


and u) to compatibility (5.1) and AMB (5-3) yields a first order and a 

second order partial differential equation for $ and u. We call this 

<•* « 

the Second boundary value problem. 

The First boundary value problem Is appealing to the practitioner 
because of the clear physical significance of the principal variable, 
the velocity field, and the simplicity of the boundary conditions. The 
Second boundary value problem is unattractive because the stress 
functions lack physical appeal, the boundary conditions are complicated, 
and because It involves two coupled equations in two unknowns. In 
either case, construction of closed form solutions is virtually 
Impossible. It is this intractability which leads us to search for 
and study methods by which the solutions of these two boundary value 
problems may be approximated. 

Two approximation methods for boundary value problems of solid 
mechanics dominate the current literature, namely, the method of 
finite differences and the finite element method. The finite dif- 
ference method is, for the most part, reserved for dynamic (i.e. wave 
propagation) problems. Only rarely is it used for analysis of quasi- 
static deformation problems. In using the finite difference method, 
one attacks the partial differential equations as they stand. However 
the finite element method requires that one first generalize the 
equations (5.1) through (5*9). In the engineering literature this 
generalization Is accomplished by finding a variational problem which 
is 'equivalent 1 to the original problem. Perhaps nowhere can a more 
lucid exposition of the fundamental variational principles of solid 
mechanics be found than in the monograph of Washizu [36]. His modus 
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operandi is adopted here. 

The essence of Washizu's method is to use Lagrange multipliers 
to 'enforce' the equations (5-1) through (5*9) of the general boundary 
value problem. These multipliers turn out to be the displacement, 
stress, traction, etc. After discussing two 'virtual work' principles, 
he supposes the existence of a potential for the stress (the strain 
energy), and shows that all the variational principles of classical 
linear elasticity may be systematically derived from the "generalized 
potential energy principle" now known as the 'Hu-Washizu' principle. 

The principles discussed in this work are named 'virtual wcik,' 
'potential energy,' ' Hel 1 inger-Reissner, ' etc., in analogy to the 
principles of linear e lastostat i cs . 

We begin with the generalization of the linear momentum balance 
(5.2). Let us momentarily regard 6v as a Lagrange multiplier. Then a 
stress rate t and a traction rate T satisfy LMB (5.2) and the stress 
principle if 

J [- V • t - pb] • 6vdV - J (T - n • t)<5ydS ■ 0 (5. 13) 

for arbitrary 6v. In (5*13), as in (5.2), t apparently must be dif- 
ferentiable. but 6v need not even be continuous. Now, by formal 
Integration by parts,* (5.13) may be transformed to 


throughout this chapter, integration by parts is formal . 
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(5.1*0 


/ [t : VSv - pb • 5v]dV - f T • <5vdS ■ 0 

V ~ s 

to be satisfied for arbitrary 6y. In (5-1*0 apparently t need not be 

• 

continuous, but 5v must be once differentiable. Any stress rate t 
admissible to the classical form of LMB (5.2) is also admitted by 
(5.1M, but the converse is not true.* Therefore we call (5.1*0 the 

m 

'generalized' linear momentum balance, and write LMB (5-1*0. Any t 
and T which satisfy the stress principle and LMB (5-2) necessarily 
satisfy (5-1*0. 

Consider the modified general boundary value problem composed by 
replacing LMB (5.2) by generalized LMB (5.1*0 and eliminating compat- 
ibility through use of (5-10): 

e * j(Vv T + Vv) ; hi - i(Vv T - Vv) ; 

f [t : V6v - pb »5v]dV - f T • 6vdS ■ 0 

- s 

for arbitrary 5y ; 

t - G(e,u),x) ; T t " I t on S o ; 
v * y on S y . (5-16) 

Every solution of the general boundary value problem (5.1) through 

*in particular, (5.1**) admits piecewise continuous stress rates, 

whereas (5-2) does not. 
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(5.9) is a solution of (5-15) and (5.16). If we admit only y-y on 
S y , 6v - 0 on S^, and T on , then one easily reduces (5.15) to 
a functional of the velocity field alone. This comprises the principle 
of virtual work. The virtual work principle is the most common basis 
for finite element algorithms used in engineering today. Its strong 
appeal stems from its overall simplicity. In arriving at (5.15) no 
assumption was made which restricts the from of the constitutive 
function G, except that all t*G(e,w,f) satisfies AMB (5-3). As an 
alternative to retaining VBC (5.16) as a subsidiary condition, we may 
'enforce' it by introduction of the Lagrange multiplier 5J or. S^. 

We replace (5-16) by 

/ 5T • (v - v)dS » 0 for arbitrary 6T . (5.17) 

S _t 

v 

If we admit only J ■ J t on , then one easily reduces (5.15) through 
(5.17) to a functional of only v and T (on S y ) . We call this the 
Second virtual work principle, it is not used as a basis for finite 
element algorithms in engineering because the velocity boundary con- 
dition is so easily enforced; that is to say, introduction of the 
extra variable is a greater effort than a priori satisfaction of 
VBC and 6v * 0 on S y . 

By a procedure similar to the one above we may obtain the 
generalization of the compatibility equation. Let us momentarily 
regard the stress function 4> and the surface tangents 6s^ as 
Lagrange multipliers. We write (Ss^e.) as (nx<54>). If e and u (in V), 

- -j ” 
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and y (on S) satisfy compatibility (5.1), then they necessarily satisfy 


/ [V x (- e +oj)] : 6 $dV + f (n x 5 $) ; (- e+w+7v)dS - 0 

"y » * * "j “ ~ ~ ~ — 

( 5 . 18 ) 

for arbitrary <54f. Just as for LHB , we formally integrate by parts in 
order to relax the smoothness requirement on e and uj, obtaining 

/ [-e +u] : (7xj*')d V + f (n x $$) : (Vy)d$ » 0 (5.19) 

V ~ $ 

for arbitrary <5$. For the same reasons that we called (5.1^) a 
generalization of LMB, we call (5 - 19) a generalization of compatibility. 

In this form it is easily seen that only the in-surface components of 
(Vv) influence the value Of the functional. If E, u), ancf y are found 
which satisfy (5.19) identically for all 5$, It Ts incorrect to con- 
clude that (e-aj)»Vy anywhere in V, or even on S; only the in-surface 
components of (Vv) will agree with (e - uj On S. Unfortunately this 
fact is not brought out in the literature, and is obscured by the 
conventional form of (5*19), which we now give. Using the formula 
(integration by parts) 

/ (n x 69) : VydS » f n • (V x 69) • vdS (5.20) 

S ' S " * 

and identifying (V x 54>) as 5t (in agreement with 5.11), equation 
( 5 - 19 ) becomes 
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(5.21) 


J [- e + to] : Std V + f n • St • vdS * 0 . 

v - $ " "" " 

It Is clear from the present derivation that St (In 5.21) is subject to 
no constraint except 

St « V x 64> . 

Consider the modified general boundary value problem obtained by 
replacing compatibility (5-1) by generalized compatibility (5-21), 
and using (5.11) to eliminate LMB: 


(5.22) 


(5.23) 

Every solution of the general boundary value problem (5.1) through 


f [- e +oi] : StdV + J n • St * vdS * 0 

v ~ ' ~ s' ~ ■ 

for all St * V x 6$ ; 
t * Vx $ + t* 3 ; n • t ■ J t ; to + to"^ ■ 0 ; 


(e + to) • x + t - t T - T'(e-to) - 0 ; 
e ■ G 1 (r,x) ; r - i(t + t • u - to • t + t T ) ; 


v * v on S 


It ’ *t on S a * 
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(5.9) Is a solution of ( 5 . 22 ) and (5-23). If we admit to (5.22) only 
v* v on S , n • t ■ T on 5_, n • <$t » 0 on 5 . and (t,ui) which satisfy 
AMB, then we can reduce (5*22) and (5.23) to a single functional. This 
functional, when we are able to construct such t and u>, comprises the 
principle of complementary virtual work. Except for pathological 
cases, construction of those t and oj is Impracticable. The problems 
associated with the use of the complementary virtual work principle in 
its 'pure 1 form (5.22 and 5-23) may be avoided if AMB (5-3) and TBC 
( 5 . 7 ) are incorporated into the principle by introduction of the 
Lagrange multipliers 6 u) and Sv: 

/ [ (e + u) • T + t] : dV « 0 (5.24) 

- - - - ~ 

for all 610 : 6 u> + » 0 ; 

-mm ‘4* T*# 

/ [n • t - T. ] • 5y dS » 0 (5.25) 

* - 1 * "t “ 

s a 

for arbitrary <5v (on 5^) . 

We call (5.22) and (5.23) thus modified the Second complementary 
virtual work principle . 

For easy reference we collect the equations of the Second 
complementary virtual work principle below: 
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/ {- e + ui] : 6tdV + f n • 6t • vdS * 0 
- ~ - J - - - 

V S 

for all 5t : 6t * V x 5<t> ; 

m* m* 


t«Vx$+t^; u> + oj"^ • 0 ; 

•'* ■v «« V ^ 


J [(e + t»j) • t + t] : 6udV = 0 
V ' ~ " ” 

for all 5u: 5w + 5w T = 0 ; 


e • G 1 (r,x) ; r“i(t+T*u>-(jj*T + t T ) 


S (o * t _ I* ) * 5yds * 0 

f L 


for all 6y on S 0 ; 


v " v on S v . 



( 5 . 26 ) 


It is a simple exercise to reduce (5.26) to three functionals 
involving t and u in V and v on S^. The second complementary virtual 
work principle used in this manner is the basis for the finite element 
algorithm presented in this work. In deriving (5.26) no assumption 
has been made which restricts the form of G‘, except that e-G'(r,T) 
must be symmetric. 

The most important feature of the virtual work principles is 
that they admit functions v, e, u), and t less smooth than did their 
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partial differential equation counterparts. A second property worthy 
of special note is that they were derived without restricting the forms 
of the constitutive functions G and G 1 . The various energy principles 

A* m0 

we are about to derive require such restrictions. Specifically, we 
require that a stress-rate potential W for r exist such that 

r * 3 W(c,t)/3e . (5-27) 

Moreover, we consider only cases in which the nominal traction rate 

T is prescribed on S a .* for solids of the type (A. 20) such a 

potential exists if and only if V (and hence W) is symmetric. We 

ax 

remark that the constitutive equations for solids common In the 
engineering literature all satisfy this symmetry.** 

If the conditions above are satisfied, then the first energy 
principle is found by introduction of the potential (5.27) for the 

stress rate to the virtual work principle (5 - 17) and (5.18). But 

■ 

since that principle involves the stress rate t, we need a potential 
U for t such that 

t - 3 U(L,t)/3L T . 

Recalling the definition of t and r, we get 

*this is actually over-restrict i ve; see Hill [37], and 
references to Sewell, therein. 

**including all the materials in examples accompanying this work. 
Some exceptions are noted by McMeeking and Rice 13*0 ■ 
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(5.28) 


i(t + t T ) * r - i (t * oi - oj*t) 

** •» •• ■v «# ^ 

and recalling AMB (5.3) we get 

i(t - t T ) - -il(e + w) • t - t • (e- u)] . (5.29) 

Putting 6L T » <5e - <$u), we get from (5.28) and (5.29) 

• T • 

t : 6L ■ r : 5e - (6e : t • w + e : t • 6w) 

- $T : (to • 5u> + • w) . ( 5 . 30 ) 

It is easily seen by inspection of (5*30) that t : 5L^" is a perfect 
differential if and only i f r : 6e is a perfect differential. When 

r : Se * (3w/3e) : $£ 

then 

t : 6L T - (3U/3L 1 ) : 6l T , 

where W is defined by (t is suppressed) 

W(E) » 1 e : W : e + e : I (5. 3D 

and U is defined through W as 
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(5.32) 


U(e,w) ■ W(e) - e : t • u - i t : (u • u) . 

•>» +0 % -M IV «v 

Now we Introduce the potential (5-32) for t to the Second 
virtual work principle (5.15" 5.17) to obtain 


e - i(Vv T + Vy) ; u - i(Vy T - Vy) ; 


(5.33) 


STTfv.E.u, J ) - 0 ; 


(5-3*0 


■rr(y,£,w,T ) « f [U(e,u) - pb»v]dV-f T • vdS 


- / f • (v - y) 


dS 


There are two ways to deal with the subsidiary conditions (5-33). It 
is an easy exercise to reduce (5-33) and (5-3*0 to a functional of 
the velocity field and traction rate. Alternatively we may ’enforce 1 
(5-33) by use of Lagrange multipliers. 

The first course of action leads us to the principle of 
stationary potential energy 


<5tt (v,T ) - 0 

p - -t 


(5.35) 


tt (v,T ) - ir(v,}(Vv^ + Vv),i(Vv T - Vv) ,T ) 
P " " t ~ ” " -**t 


Any solution of the general boundary value problem is a solution of 
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(5.35). In practice one usually satisfies the velocity boundary 
condition a priori, thereby eliminating T from (5-35). 

The alternative course of action leads us to a Hu-Waehizu 
energy -principle. Let t be regarded momentarily as a Lagrange 
multiplier. Then (5.33) and (5.^) may be set as 

6TT HW (y,e,w,j t ,t) = 0 ; (5.36) 

“"hv/- ’JE’-'-t’ £) " ^(y.e.u), J t ) + f t : (Vy - (e-w))dV 

e - e T - 0 ; u + » 0 . 

It is possible to enforce both of the conditions (5-33) with the single 
multiplier so long as only symmetric e and skew symmetric w are 
admitted to There are no other subsidiary conditions. Any 

solution of the general boundary problem (5.1 through 5.9) necessarily 
satisfies (5-36). We write out the stationary conditions for future 
references: j_ 




LMB: 


f [t : V6v - pb • 6v]dV -/ T • 5vdS - f T • 6vdS - 0 (5-37) 

V ~ ' ' S„ S " 

a v 
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CONSTITUTIVE EQUATION: 


f [ OW/3e) - i(t + T»o)-u*T + t T )J : 6edV -0 (5.33) 

V 


AMB: 


f [i( (e + co) • t + t - t T - t • (e - w) ) : 6o)]dV * 0 (5-39) 

V 


VBC : 


J 5T • (v - v)dS ■ 0 (5.^0) 

S v" 

COMPATIBILITY: 

/ [7v - (e- u)] : 6tdV - 0 . (5.*»1) 


Notice that the stationary condition for 5v (5.37) is the generalized 
LMB (5.1*j\ 

The main detraction of the Hu-Washizu principle is the large 
number of unknowns (five). Although a finite element algorithm could 
be based upon this principle, it is unlikely that it could be made 
very efficient, and thus would be of diminished practical interest. 

We note that solely by rearrangement of terms, (5.26) may be 
rendered in the form [15] 
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' Fr HW (v,e,a>,T t: ,t) 


( 5 . 42 ) 


f [W(e) - i(t + T»co-w»T + t T ):e] - lx: wu + t : u 

m* — ^ <V ** »* rn» «. J 

f |t : Vv - pb*vjdV - J • vdS - J J t • vdS 


dV 


+ f T • v dS . 
S _t 


I f the constitutive equation (5-5) is used to eliminate z as a variable 
from (5.*»2) , then, defining 


-R(t,o)) » W(e(t ,uj) ) - i(t + T*u-u*T + t T ) : e(t ,u>) , 


we obtain a Hellinger-Reissner variational principle: 




(5. A3) 


xr HR (v,u),j t , t) “ J |-R(t,u>) - i x: (u • u) -t- t : wjdV 


+ f T • vdS 

S _t 
v 

+ f [t : 7v “ pb • v]dV - f T • vdS - f T • vdS ; 
V ~ ■ s ’ 1 " s 


w + ui * 0 . 
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Any solution of the general boundary value problem is also a solution 

of (5 • ^3 ) - The stationary conditions are LMB (5.37). AMB (5.39), VBC 

(5*40), and compatibility (5-41). In the stationary conditions e only 

appears as a function of t and u. 

If we admit to tt (5.43) only stress rates t and traction rates 
HR — 

T which satisfy the generalized LMB (5.14) for assigned T , we obtain 
the complementary energy principle of Atluri [15]: 

6tt (w,t) « 0 ; (5*44) 

c — 

tt (ui,t) - f }-R(t,u)) - (w ui ) + t:tojdV + f n • t • v dS 

C ~ V ~ ~ S 

v 

• • b T 

t * V x $ + t ; u + w *0; 

•St *V “w 

n • t * T on . (5.45) 

Any solution of the general boundary value problem Is a solution of 
(5.44) and (5.45). The stationary conditions are AMB (5. 39), VBC 
(5.40), and compatibility (5-41). In the stationary conditions £ only 
appears as a function of t and u). It is not a simple matter to reduce 
(5.44) and (5.45) to a single functional of t and u since this would 
require construction of stress functions satisfying 

n • (V * $+ t b ) » t. on S„ . 

* % * t V 
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Therefore (5.44) and (5.45) are not suitable as a basis for a finite 
element algorithm. 

The problems associated with the use of the complementary energy 
principle In Its 'pure' form (above) may be avoided if the traction 
boundary condition is incorporated into the principle. Let v be a 
Lagrange multiplier on S^. Then we replace (5.44) and (5.45) by 

5TT*(v,ta),t) * 0 ; (5- 46) 

■n' (y.cj.t) ■ tt (u,t) + f (n • t - 7 ) • vdS 

C " ^ C -w » g * — “ t ” 

G 

• • b T 

t * Vx$ + t ; u + u * 0 . 

Any solution of the general boundary value problem is a solution of 

(5.46) . The stationary conditions are the same as those of the 'pure' 
complementary energy principle, except that the traction boundary con- 
dition follows from dTT^/Sy - 0. It is a simple exercise to reduce 

(5.46) to a single functional of v, o>, and t. The complementary energy 
principle thus modified serves as the basis for the finite element 
algorithm presented in this work. 
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CHAPTER V! 


UNIQUENESS OF SOLUTIONS OF BOUNDARY 
VALUE PROBLEMS 

Int roduction 

Once a solution of a boundary value problem has been found, one 
may Inquire as to the uniqueness of that solution. In the event that 
multiple solutions exist, a stability criterion is required to dis- 
tinguish stable solutions from those which are unstable. 

Hill [38] has published a number of papers on the subject of 
uniqueness. For materials of the type (4.20) through (4.22), his 
criterion is necessary and sufficient for uniqueness. He has also 
given a practical method for testing the uniqueness of deformations of 
hypoe 1 as t i c/plas 1 1 c bodies (which are nonlinear in a certain way). For 
extended treatments of the theory, see Hill [38]. 

In this work we establish (quas i stat i c) stability of a configura- 
tion of a body by Inspection of the 'restoring force' that results when 
the configuration Is slightly perturbed. If the forces arising from an 
admissible* perturbation tend to attenuate the disturbance, then we say 
the given configuration is stable. The principal shortcoming of this 
stability criterion is apparent when uniqueness is lost; It provides no 
means of distinguishing stable solutions from unstable solutions. The 


* 


not violating kinematic constraints. 
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criterion tests configurations of bodies, not solutions of the boundary 
value problem. This shortcoming is probably unavoidable when one is 
using a purely mechanical material model such as (k.20), or the 
hypoelastic/plastic model. 

We construct two uniqueness criteria, one from the virtual work 
principle, and the other from the Second complementary virtual work 
principle. These are equivalent to the criteria proposed by Hill [38]. 

A Uniqueness Criterion Based on the 
Virtual Work Principle 

Let (y\ t') and (y 2 , £ 2 , u 2 , t 2 ) be two solutions of 

the general boundary value problem. The same traction boundary con- 
dition and velocity boundary condition are satisfied by each of the 
two solutions above. By the virtual work principle, any solution of 
the general boundary value problem is also a solution of ( 5 . 15 ) and 
(5.16), so (v' , , w* , t') each satisfy (5.15) and (5.16). Therefore 

the difference of the two solutions 

(Av,A£,Aco,At) ■ (v 2 - v 1 ,e 2 - e 1 ,w 2 - w 1 ,t 2 - t 1 ) (6.1) 

necessarily satisfies* 

Ae - i(VAy T + VAv) ; Aw - i(VAy T - VAv) ; 


the first of (6.2) applies when nominal tractions are prescribed 
on S^; the second applies when true tractions are prescribed. 
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or 


( 6 . 2 ) 


/ IAt : 7<5v - pAv • Vb • 6v]dV ■ 0 ; 

V ~ ” 

fy[ At : V5y - pAv • Vb • 6y]dV - J [ (A J - n * Ae • n) J • 6y]dS * 0 ; 

a 

At ■ G(e 2 ,w 2 ,T) - Gte 1 , a) 1 ,t) 

v * 0 on S : 

- v ’ 

A sufficient condition for uniqueness is therefore that (6.2) have no 

solution among all pairs {v ,v } taking the prescribed value y on 5 y . 

Generally we are unable to reduce the expression above to a functional 

1 1 2 

of Av alone since the constitutive equation involves £ , W , E , and 

" *v 

2 

u distinctly. However, if G has the distributive property (modulo 
G(0,0 ,t)) 

(6.3) 

then it is an easy exercise to reduce (6.2) to a functional of the 
velocity field alone, and that functional is independent of the 
particular boundary values v. This is important In practice because it 
means that searching for solutions of (6.2) is precisely equivalent to 
searching for solutions of virtual work (5.15 and 5.16) with homogeneous 
boundary conditions and no relaxation. 

The condition (6.3) is satisfied by hypoelastic and hypoelastic/ 
viscoplastic materials, but not by hypoelast ic/plast ic materials. 


G(Ac,A(jl),t) - G(eV.T) - G(e .u'.t) + G(0,0,t) 


5v ■ 0 on S 
- v 
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Hill [38] has shown that if the 'linear comparison solid' associated 
with a particular plastic body is used in (6.2), then uniqueness of 
solutions for the linear comparison solid is sufficient for uniqueness 
of solutions for the plastic body. 

A Uniqueness Criterion Based on the Second 
Complementary Virtual Work Principle 
In this section we use the same difference notation as in the 
previous section. Any two solutions of the general boundary value 
problem necessarily satisfy the Second complementary virtual work 
principle. Therefore, their difference satisfies 


f [-Ae + Aw] : 6tdV + f n ♦ <5t • Avd$ - 0 

II m* *» » r 


for all <5t * V x 5$ ; 


At ■ V x A$ ;* Aw + Aw * 0 ; 


j [ (Ae + Aw ) • T + t ) : 6wdV ■ 0 
V ~ ~ 


for all 6w: 6w+5w *0 ; 


- ' ' *2 . -v - ' / • 1 

u \ r , .; - G (r ,T) ; 


r ■ r(t+T •u-w*t + t T ) ; 


(6.M 


for convenience, we have assumed that Ay • Vb * 0. 
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In this case searching for solutions of (6.4) is precisely equivalent to 
searching for solutions of complementary virtual work (5.26) with 
homogeneous boundary conditions and no relaxation. Of course, the 
same materials satisfy the condition (6.5) which satisfied the con- 
dition (6.3). In the case of hypoelastic/plastic materials, (6.5) is 
not satisfied, but we note that the constitutive function of the 
associated linear comparison solid for a particular hypoelastic/plastic 
solid does satisfy (6.5). In view of the fact that uniqueness for the 
linear comparison solid Is sufficient for uniqueness for the 
hypoelastic/plasti c solid (and no restrictions are placed on the 
method by which we establish uniqueness for the linear comparison 
solid), we may use the constitutive equation of the linear comparison 
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solid in (6.4) to establish uniqueness for the hypoel as t i c/pl ast i c 
solid. A proof of this sufficiency would be difficult if one started 
from (6.4), since plasticity theory is not formulated In terms of the 
stress rate r. In the bifurcation study accompanying this work, the 
criterion (6.4) is used in conjunction with the constitutive equation 
of the linear comparison solid. 
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CHAPTER VII 


A FINITE ELEMENT ALGORITHM 
Introduction 

The finite element method had its beginnings in structural 
analysis, but spread quickly and with great success to other areas of 
applied science. The underlying mathematical theory of finite elements 
resulted from the study of the method as applied by engineers, and is 
still being developed at this date. Extended introductions to the 
finite element method, both from the practitioner's point of view [39], 
[LO], and from the mathematician's point of view [LI], [L2], are widely 
available, and therefore omitted from this work. For the discussion 
that follows, the finite element method may be regarded as a generaliza- 
tion of the approximate methods based on energy principles of linear 
elastostat i cs . * 

The finite element algorithm described in this chapter is based 
upon the Second complementary virtual work principle (5.26). The 
advantage of starting from the 'work' principle, as opposed to starting 
from the 'energy 1 principle (5-L6), is not only greater generality with 
regard to constitutive equations and boundary conditions, but greater 
clarity, since each of the equations of (5.26) corresponds to an 
equation of the general boundary value problem (5.1) through (5-9). 

* 

see Washizu [36], especially sections 1.5 and 1.7 therein. 
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When the constitutive equation and boundary conditions are such that a 
complementary energy principle exists, the distinction between the 
'work' and 'energy' formulations vanishes entirely. We do assume that 
the constitutive equation may be set in the form (5.5), and the traction 
boundary conditions are dealt with as if nominal traction rates were 
prescribed, but generalization to treat other materials and types of 
boundary conditions should present no difficulty for the reader 
familiar with finite element methods. 


A Finite Element Algorithm 

According to the Second complementary virtual work principle 
(5.26), the stress rate t, spin u, and boundary velocity y of a body 
are solutions of the following boundary value problem: 


y j[(W 1 : (r - E) + u>) • T + t J : 6uj 


dV « 0 


J * j |-W ^ : (r -E) +uj:6t|dV + J" n • fit • v 


vdS “ 0 


(7.D 


for all <5u): 5u+6u ■ 0 ; 

6 t: 6 t ■ V x 5$ ; 


J (n • t - J t ) • 5y 


6vdS ■ 0 


for all 5v: / <5v • 5vdS + 0 

■'s 


(7.2) 
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with the subsidiary conditions: 


• *o *b *o _ *k • 

t ■ t + t i t * V M ; V • t D ■ ~pb : 


• • • ▼ ▼ 

r»i(t+T*w-u»T + t ) ; a> + u> * 0 ; 


v » v on S 


In formulation of a finite element algorithm we regard the body 

as an assembly of sub-bodies called elements. We wish to represent the 

functions t and id independently on each element. The stress rate t will 

be represented as indicated in the subsidiary conditions above on 

•o 

each element, but between elements t will generally be discontinuous. 
Such t is still admissible to the complementary virtual work principle 
providing it satisfies the generalized LMB (5.1*0.* Indicating by 
NELM the number of elements into which the body has been partitioned, 
we write LHB (5. 1*0 as 



T • SvdS • 0 
-t 


for arbitrary 5v (continuous across interelement boundaries) 


Integration by parts yields 


note that the original form of LHB C5.2) does not admit such t. 
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where and S N are the element domain and boundary, respectively. 

• JL 

Finally, setting T to T on 5 and 6v to zero on 5 , we obtain 
3 -t -t a - v 


!£k M f r . r i 

/ 4 j I (n • t • 6v)dS - / J t 

N= 1 c .. f c e ^ ( c i 


• SvdS > = 0 


VW 


(S„ ns ^) 
n a 


(7.3) 


We replace the traction boundary condition (7.2) by this last equation 
(7-3). which is easily seen to be no more than a statement of traction 
reciprocity. In the context of finite element terminology, it is 
called 'interelement traction reciprocity.' One should take special 
note that in (7.3) Sy is required to be single valued on the inter- 
element boundaries, since 6v is required to possess a (generalized) 

gradient everywhere in V. 

• • 

Since t, w, St, and Su are now independent on each element, 
(7-1) must be satisfied independently on each element. Thus, we 
replace the boundary value problem (7-1) and (7.2) by 



{ (W 1 : ( r - Z) + w) • T + t ] : 6o) } dV » 0 


f-W * : ( r - I) + w] : St JdV + f n • St 

= - - - -J v - 

for all 5w: Sw + 5uJ * 0 ; 

•%» *v ■** 

for a 1 1 5t : St ■ V x 6$ 


ydS 



(7.M 


79 



NELM / 

2 f (D*t-5y)dS - J f t 


6vdS > - 0 


(7.5) 


N ” S N' tS N n V 


(S N n V 


with the subsidiary conditions 


• *0 # b *0 • 

t » t + t ; t - V M ; V • t D - - pb 


T • • • T 

u) + (0 «0 r ■ J (t + T # u)-w # T + t ) ; 


on each element, and globally 


v = v on S y ; 5v - 0 on $ y ; 


y and 5v single valued along interelement boundaries. 


To construct a finite element algorithm based on (7.**) and (7-5) 

• • 

we must be able to find representations for v, fiv, u), t, fiw, fit, and 

* m n* »v -w 

r which satisfy all of the subsidiary conditions explicity. In the 
next few paragraphs such representations are discussed. 

Let us represent the velocity on the surface of the N th element 
by 


i ( *> *2 

i-1 


<i 


- , _ _ 
q on S U DS 
N v 

q elsewhere 


(7.6) 
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where the N. are (vector valued) Isoparametric shape functions* and q 1 
are nodal velocity parameters. The shape functions are NQ in number. 
Those parameters determined by the velocity boundary condition are 
distinguished from the undetermined parameters by an overbar. The 
variation 6v may be found from (7-6) as** 

NQ 

<5y(x) - Y! N.WSq^ • (7-7) 

i«l 

Note that v * V and 6v » 0 on according to (7.6) and (7. 7). Further- 
more, when using the isoparametric shape functions, we can assure that 
v and <5v are single valued along any interelement boundary simply by 
'connecting* the nodes of the elements adjacent to that boundary. The 
easiest way to deal with this connectivity in practice is to index the 
velocity parameters node by node, for the whole body, instead of having 
indices which are independent on each element. The relation between 
the parameters with element-level indices (q^) and those with global 
indices (Q 1 ) is formalized by the introduction of an 'assembly matrix' 
[A u ] for each element such that 


*see Ergatoudis et al. [A3]. 

It Is unnecessary for v and Sy to be so related when starting 
from a 'work' principle; similar statements may be made for w and 6 uj, 
t and 6t. We relate the total quantities to their variations so that 
'work' and 'energy' formulations will coincide when an energy formula- 
tion exists. 
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{q N > - [a n H5) ; 

and 

- [A n H5q}. 

It turns out to be more natural to deal with the velocity boundary con- 
dition by constraining Q than by constraining q^, so the q^ and the q N 
(see 7.6) never need to be distinguished explicitly in practice. The 
total number of parameters Q* is called the 'number of degrees of free- 
dom' of the finite element mesh, abbreviated NDOF. 

We represent the spin u and the stress function <f> in the N th 
element's Interior by 

NW 

w(x) QW. (x)a^ (7.9) 

i = l 


Q on S. 


Q elsewhere 


(7-8) 


«(*) *2%i (i)S N ' 

i*I 

QW. and 0. being tensor valued shape functions. The actual functions 
used in the particular examples accompanying this work are detailed in 
Appendix C. The shape functions for the spin are chosen so that each 
satisfies , -- r --. ^ ... 

QW. + qwT * 0 . (7.10) 
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The stress rate Is represented in each element by 


NT 

t(*) “ £ <£; ( x) 6 |!, + t b <>0 (7.11) 

i -1 

where 


qj. (x) « V x <J>. (x) ( 7 . 12 ) 

for each parameter 8 ^. From (7.9) and (7.11) we obtain representations 
for 5w and St as 


NW 

6w(x) - 53 ^i ( ? )6a N 
i -1 

NT 

s£W * £ Q Ji { ? )66 n 

i = l 

• • • T\ 

Finally we form r = i(t + T*u)-u) , T + t ) as 

mt «» ** — * 

NT # NW 

r(x) «53 i(< £i +( ^I )e N + £ * ( I * 9^i ‘ Q ^i -lK * (7-lM 

i=*l 1 = 1 


(7.13) 


• • • 

The functions v, 6 v, w, t, 5u, St, and r, when represented in the forms 
( 7 - 6 ), (7.7), ( 7 - 8 ), (7.9), (7.H), (7.13), and (7.1*0, satisfy all of 
the subsidiary conditions previously mentioned, and are therefore 

admissible to the boundary value problem as stated by ( 7-*0 and ( 7 - 5 ). 

• • • 

Putting the representations for y, ov, u, t, Sw, St, and r into 
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the functionals (7-M and (7-5) and carrying out the assigned 
integrations yields the following finite element equations (in which 
the index * N ' Is suppressed on the stress and spin parameters): 


(6a} T | 

: -[h 11 

h' 2 j|“| * ( p a - b ) 

+ {p 01 * 1 }}- 0 

(7.15) 

{66) T 

! -[H 21 

h 22 ]|“s * ( p B - b ) 

+ {p 8,I > 

(7.16) 


+ [G]{q N }j= 0 


“tLn, _ , 

ZK>“> B - 1- W - ° • (7 -' 7) 

N=1 

f 

Henceforth we refer to (7-15) as AMB, to (7-16) as compatibility, and 
to (7-17) as TBC. The individual matrices are defined below: 


h!! = f {(t • QW.) : D : (t • QW.) + t : (QW. • 0W.)}dV 
•J *\i ~ •» ~ '•J 1 J 

N 

(7.18) 

h! 2 = / {(t • QW.) : D : (QT.) - QW. : QT.}dV 
V N 

(7.19) 

H 2 ! « / {(QT.) : D : (fQW .) - QJ : QW }dV 
N 

(7.20) 

H? 2 - / {(OT.) : 0 : (0T.)}dV 
N 

(7.21) 

G. . - / n • (QT.) • (N )dS 

(7.22) 

1 J 1 J 

5 n 



84 




(7.23) 


F ! - J ?,•«¥,)« 


p?- b . 
1 


j(Q IV.) : [t b + (D : t b ) • x + ^ t] JdV 

(7.2k)* 

p B,b „ 


|(0T.) : (-D : t b )JdV 

(7.25) 

P?’ Z » 

i < 

i, i 

N 

(t • qw.) : D : Z JdV 

(7.26) 

i 

K 1 

N 

|(QT.) : D : zjdV 
1 1 - -} 

(7.27) 


and D Is obtained from W ' by symmet r i zat i on : 


L . . . * i(W. ! . . + W. , + W.1 1L + W . ! . ,) 

i jkl i jk 1 j ikl i j Ik j i lk 


(7.28) 


This symmetrizatlon is easily done after W is computed, and serves 

2J 

to reduce by a factor of four the number of multiplications required 
to compute the H matrices (7.18) through (7.21), and other matrices 
Involving W The integrations must be performed numerically since 
the integrands and the domain of the element change during the deforma- 
tion. in tne cXoiiip i es accompanying this work only quadrilateral 
elements were used, so symmetric Gaussian quadrature rules were used 
for the i ntearat ion. ** 

"the last term in the Integrand is a residual whose significance 
is explained in the next chapter. 

•frh 

see Tor.g and P.ossettos [39] t especially chapter six. 
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The procedure which leads one from equations (7.15) and (7.16) 
to the element stiffness matrix is virtually identical to that of Plan 
13]. We define the element 'H-matrix 1 as 


[H] 


H 21 



(7.29) 


b T 

and loads (P } and (P }, due to body force and fluidity, respectively, as 



(7.30) 


Then (7.15) and (7.16) may be collected into a single equation as 



+ P Z }. 


(7-31) 


If D is symmetric, that is, if W. ! k , » w k ] . j , then from (7-18) through 
(7.21) we easily determine that [H] is symmetric. 

If the H-matrix is not singular, then we solve^the matrix 
equation (NQ + 1 right hand sides) 


[H] [H~ 1 G 


H _ 1 P ] - 




(7-32) 


on each element. Explicit calculation of the inverse of [H] is not 
only unnecessary, but substantially Increases the time required to 
generate the element stiffness matrix. According to (7.30. The spin 
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and stress parameters on each element are given by 


[‘i-rn-W + ^P). (7.33) 

Using (7.33) to eliminate {a/S) from TBC (7.17) leads to 
NELM 

S ! { v T iw + {s v T t° G N) { H''p n ) 

N-l 

- ■ ° (? - 3M 
in which the element stiffness matrix has been identified as 

IK n ] « [0 gJ] [H H G n ] , (7.35) 

and the resultant nodal 'forces' are given by 

- [0 gJ] + (F n ) . (7.36) 

It is easily verified that the element stiffness matrix [K] is symmetric 
if [H] is, and so the symmetry of [K] ultimately depends upon the 
symmetry of the constitutive matrix W. 

To this point all of the finite element equations are indepen- 
dent on each element. The formal assembly of the global stiffness 
matrix and loads is accomplished by introduction of the assembly 
matrices (see 7.8) to (7.3*0 so that the element level velocity 
parameters may be expressed as functions of the global velocity 
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parameters. For {q } and { 6q } we write 


{q N > ■ IA n ]{Q + 0) ; {6q N ) - £A n ] {6Q} ; 

and from (7.3 1 *) thus obtain 

{SQ} T [K g Hq} - {6Q} T {P q } - {6Q} T [K G ]{Q} . (7.37) 


In equation (7*37) the global stiffness matrix [Kg] and the loads 
{ P G } are defined by 


NtLn 

£ 'v "v'v 


(7.38) 


and 


{P G > « [a/ j (f n ) - [0 ^]{h-'p n } . 


(7.39) 


The load matrix {P-} contains contributions from the prescribed 

G 

body for rate b, the relaxation Z, and the tract ion boundary j:on- 
dition T . The global stiffness matrix, as defined by (7.38), will 
be singular for rigid translations, but not for rigid spin (except .in 
the case that there is an 'axis of equilibrium' [A4]). In order to 
solve the equation (7-37) we define a modified global stiffness matrix 
[K*] and a modified load {P*} as follows: 
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1 1 Hill 1 



(7.40) 



S U if ( «l ' V or <5j - Qj) 

Kjj otherwise 


P 


* 

I 



if (Q, - Q, ) 

NELM 

- K | jQj otherwise . 
J=1 


(7.1*1) 


Then (7-37) may be replaced by 

[ K*] {Q} * {P*} . (7.1*2) 

If [K*] is not singular, then we solve (7-1*2) for (Q.) , 

{Q} - [K*f ] {P*} . (7.1*3) 

By backsubsti tutlon we obtain the velocity (on the boundary of each 
element), the spin and the stress rate on each element: 


(q N > - [A n ][k*3 1 (P*> 

(7.1*1*) 

|^J " I H " lG N nAN3[K*j” 1{ P*} + {H ’ lp N } 

(7.1*5) 

v(x)» |N(x) 1 [A n ] [K*]" 1 {P*} 

(7.1*6) 
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0 

QT(x) 


|[h' 1 G n ][A n ]IK*3‘ 1 {P*} 4 {H 



(7. *7) 


Equations (7-^*6) and (7.^7) comprise the approximate solution of the 
boundary value problem. 

The velocity v(x) is determined by (7.^6) on the element 
boundary only. On an element's interior we may construct e(x) from 
t(x) and u)(x) according to the constitutive equation (5-5). If our 
computations have produced an exact solution, then we can find the 
velocity field on the interior of an element by Integration: 


v (x ) » v (x ) + f ” (c + u) • dx , 

- -o - 

-o 

where x q is some point on the element's boundary. However the finite 
element algorithm generally will not produce exact solutions, so the 
integral expression above is of no use In defining the velocity on the 
interior of the element. Previous researchers using complementary 
work and energy based finite elements do not mention this problem, and 
it is left for the reader to assume that they found the velocity on 
the Interior of the element by interpolation of the boundary velocities 
(by use of the isoparametric shape functions). If the reader will 
recall the discussion surrounding the derivation of the generalized 
compatibility (5.21), it is evident that there is reason to doubt the 
validity of this procedure. In the absence of supporting arguments. 

It amounts to assigning the velocity on the interior of the element In 
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an arbitrary manner. 

A Jittle light is shed on this problem by the following 
heuristic argument. In preparation, we consider rules by which fields 
v, 03, and t may be judged 'admissible* or 'inadmissible* on a single 
element. Admissible velocity fields are those which satisfy 

J {v • v + (Vv) : (Vv)}dV < 00 ; (7.^8) 

V 

and the admissible velocity fields comprise an inner product space R v « 
Each velocity field in R y takes on certain values on the boundary of 
the element, and the space made up of those boundary values we denote 
R°. In a similar manner we define the space R as the space of all 

V 03 

tensors which satisfy 

f (oc : 03) dV < “ ; f (oj +oj T ) : (o3 + 03 T )dV - 0 ; (7.^9) 

V ~ ~ V ~ ~ - ~ 


and 


the space R to be the space of all tensors which satisfy 


f (t : t)dV < ® ; f (t ; V6v)dV « 0 
V " ~ V ~ ‘ 


(7.50) 


for arbitrary -Sv, 6v**0 on S. An inner product on R^ x R^ Js given by 


'“'“'r XR. 

ui L V 


f (u) : 03 1 + t : t ' )dV 


(7.5D 


where u and ware the ordered pairs of tensors 
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u » <o,t>ER x R *, 

~~ CO t 


w * <coit'>ER x R k 
•. « (*) t 


Now we define a linear map from R° x R x R to R x R by 

r V CO t (0 t 1 


w * B (u) (7.52) 

where u »<v,co,t> is any element of R v x R w x an{ ^ 

w ■ <H(e + “)*t+ t - t T - x • (e - w)] , *£+u + Vv> (7*53) 

■w <•* m* m* <¥ “ 


is an element of R x R^ determined uniquely by u. In (7-53) we have 
co t 

written £ for the constitutive function (5.5). The variational 
problem (7> i 0 may now be stated compactly as 


find u E R° x R x R. such that for all 6wER x R. 

V CO t CO t 


(6w,B(u) ) 


R X R 

CO t 


0 


(7-5*0 


(that is, 7-5*» is equivalent to 7-*0. The problem in (7*52) of 
defining w uniquely for assigned u is solved if we can show there is 
a unique y in R y such that (-e + co + Vv) belongs to R^. We note in pass- 
ing that this problem is Identical to the problem in the finite element 
algorithm of defining y on the Interior of an element. The criterion 
that (-£ +(i)+ Vv) must satisfy to be an element of R t is given by 
(7.50): 
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(7.55) 


f {(- e + w+ Vv) : V5v}dV « 0 
V ~ ~ “ 

for arbitrary 6v, 6v«0 on S. Equation (7-55), along with the boundary 
values of v (the first component of the argument of B) , and the values 
of e and oi (from the second and third components of the argument of 
B) constitute a generalized Dirichlet problem, whose solution we pre- 
sume to exist and be unique, and to depend linearly upon the argument 
of B. We conclude by remarking that the variational problem (7.5*0 
would not make sense if v were assigned in an arbitrary manner. 

The finite element counterpart of equation (7*55), including 
inhomogeneities from t and Z, is given by 

re N ){q> - nj - tc N Hq N ) ♦ {L p } . (7.56) 

where the individual matrices are defined below: 

(7.57) 

(7.58) 

(7.59) 

( 7 . 60 ) 
( 7 . 61 ) 



l? . - f VN. : [D : (QT.)]dV 
<J 

V N 


L? - / VN. : [D : (t b - I) ]dV . 
• •'u " ~ 

N 
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In equations (7*57) through (7.61) the N. are polynomial shape functions 
which vanish on S^; that Is, shape functions for 'internal' nodes, and 
the Np , and QTj are the shape functions for boundary velocity, 
spin, and stress rate which have been previously discussed. Using 
(7.1*5) to eliminate {a N /B N ) from (7-56) leads to 

{q N l - l£ M r' [uj lShh' 1 g h ] - IC N ) ia n ]{q} (7.62) 

+ [C N ] 1 1 l m l n H« p n* + * l n j| • 

Equation (7.62) expresses the velocities of the internal nodes of an 
element as a function of the velocities of the boundary nodes. In the 
sense that the problem of determining the velocities of interior nodes 
is the inverse of the 'condensation' problem encountered in ordinary 
velocity-based finite element algorithms, we might call the procedure 
above 'inverse condensation.' 

In view of the fact that use of the inverse condensation is 
potentially costly (because of the extra computation), an element for 
which Interpolation of the boundary velocity to the Interior is con- 
sistent with Inverse condensation would be highly desirable. Consider 
the following examples in which result of inverse condensation Is 
? 1 1 ustrated. 

In Lite application of a stress-based finite element algorithm 
to beam problems,* usually one can only find a piecewise linear 

*se= Murakawa et at. [ 1 3] » particularly Figure 2 therein. 
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approximation for the elastic curve, an approximation with 'corners' at 
the Interelement boundaries. The only way to Improve the approximation 
is to increase the number of elements per unit length of beam. Appli- 
cation of the inverse condensation procedure leads to (ignoring the 
length change of the beam) 

v 2 , n - -M/EI 

a familiar formula from beam theory. if the moment M is piecewise 
linear and continuous at the interelement boundaries, then oi.e may find 
a cubic spline for the elastic curve without any increase in the num- 
ber of elements. 

In the application of a stress-based finite element algorithm 
to two-dimensional problems, one typically uses four or eight noded 
quadrilateral elements. A simple example such as provided by beams 
is not available In this case, but we note that the lowest order 
polynomial function which vanishes on the boundary of a quadrilateral 
is four; being of the form 

(x 2 - 1 ) (y 2 - 1) . (7.63) 

The immediate conclusion is that the highest order complete polynomial 
which may be represented exactly on a quadrilateral with boundary nodes 
only i s' three. The shape functions for four and eight noded quadri- 
laterals involve no fourth order terms, so one could expect to gain 
nothing in accuracy by adding an extra degree of freedom for the shape 
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(7.63). On the other hand, the shape functions for a 12 (or higher) 
noded quadrilateral do contain fourth order terms, so if terms of the 
form (7.63) are Ignored, It amounts to assigning the velocity on the 
Interior of the element in an arbitrary manner. Thus, the 'inverse 
condensation' is not necessary for the 'low order' four and eight 
noded elements, but should be used if 'higher order' elements are used. 
Similar arguments may be given for triangular elements and three 
dimensional elements. In the examples accompanying this work, only 
four and eight noded elements were used. 

Numerical Stability Criteria 

The finite element algorithm just described does not necessarily 
yield an approximation to the solution of the boundary value problem. 

In reviewing the development, we surmise that the algorithm may be 
carried through to obtain the approximate solution (7>**6) and (7.**7) if 

W is nonsingular In V (see 5.5) (7.6*0 

Is nonsingular on each element (see 7.32) (7.65) 

iGuHq..} *{0} only for rigid translations (7.66) 

[K*] is nonsingular (see 7. **2) . (7.67) 

The first of these is satisfied by material models found in the 
engineering literature except for Isolated states of stress. The last 
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of these Is equivalent to the uniqueness criterion based on the Second 
complementary virtual work principle, presented in Chapter VI. 

Satisfaction of the second and third conditions above depends 
strongly upon the particular functions Nj , QW. , and QTj of the 
representations (7*6), (7-9), and (7.11). In this section we discuss 
criteria whose fulfillment is necessary for the satisfaction of (7-65) 
and (7-66). Such criteria are called ’numerical stability criteria.' 

An analogue of the condition (7-66) arises in stress-based 
finite element algorithms in linear elastostatics . The analysis of 
Tong and Pian [k5], with a minor modification, applies in the present 
case. The rank of the matrix [G] is usually 

mi n (NT, NQ - T) (7-68) 

where NT is the number of stress rate parameters, NQ the number of 
velocity parameters, and T the number of translational degrees of 
freedom of an element. It is well known that if NT<NQ-T, then 
'kinematic modes' (deformations to which the element offers no 
resistance) will occur. The rank condition which is necessary (and 
usually sufficient) for the satisfaction of (7.66) is 

NT ^ NQ - T . (7.69) 

In the examples accompanying this work, the number of stress rate 
parameters always equalled or exceeded the number of velocity parameters, 



and no kinematic mode was encountered. 


A second type of kinematic mode can occur when the finite 
element algorithm is set up In a coordinate system which has 
singularities, such as cyllnderlcal coordinates. This second type of 
mode occurs when 

NjdS ■ 0 everywhere on (7.70) 

for some particular velocity shape function. An example is provided 
by an eight noded quadralateral element In cylinderical coordinates, 
with one edge along the z axis (r*0). The shape function for the 
middle node on that edge vanishes everywhere on S except on the edge 
where dS is zero. Thus, the column in [G] corresponding to that node 
consists entirely of zeros. In such a case the kinematic mode may be 
avoided by eliminating the offending node entirely, or its value may 
be found by the inverse condensation procedure. 

The condition (7-65) turns out to be the source of most of the 
difficulty of using the present stress-based finite element algorithm. 
Even when [H] Is not singular, it may be so ill conditioned that an 
accurate solution of the matrix equation (7.32) can only be found by 
scaling.* In any case, the problem may be overcome by replacing QW. 
and QT. in a trial and error process until nonsingular [H] is found. 


that is, adjusting the magnitude of the stress and spin 
functions to improve the condition of [H]. 
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Murakawa [10] gave the necessary (but not sufficient) rank condition 
that the number of stress functions exceed the number of spin functions. 

After a number of trials, it became apparent that when [H] was 
singular, the spurious eigenmode consisted of a pure (but inhomogeneous) 
spin. Moreover, if a combination of functions [QW] and [QT] was found 
to be acceptable in the stress-free state, it remained so as the 
deformation progressed. Setting the initial stress to zero, the cri- 
terion sufficient for no spin mode to occur follows as 

[H^]{6a) / 0 ; [H^ 1 ].. «/ (QT. :QW.)dV . (7-7D 

O O Ij J — * -J 

N 

A criterion similar to (7-71), but for a finite element algorithm of 
elastic membrane theory, was given by de Veubeke and Millard [6], but 
their conclusions differ slightly from our own. A necessary condition 
for the satisfaction of (7-71) Is that 

NT* > NW (7.72) 

where NT* is the number of stress shape functions QT. whose skew parts 

do not vanish, and NW is the number of spin functions QW. . In [6], 

^ 1 

the authors suggest that the polynomial degree of the spin field m be 
related to the polynomial degree of the stress field n as 

m ■ n - 1 , (7. 73) 
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In numerical 


calling the case m^n the "classical equilibrium model." 
experience with the present finite element algorithm we find that if m 
is less than n, then the angular momentum balance is not satisfied 
with reasonable (pointwlse) accuracy. The degree of the spin field 
and the degree of the stress rate field were always the same in the 
examples accompanying this work. This amounts to the condition 

[H^] T {5B} - 0 

only for { 60 } which produce symmetric stress rate fields. 


(7.7*0 
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CHAPTER VIII 


INTEGRATION OF THE MOTION OF THE BODY 
Introduction 

Consider a body In a configuration C(t). We suppose that a 
particular finite element subdivision of the body has been defined, 
and that the stress x is known at each quadrature point on each 
element of the body. Then for assigned body force rate b, nominal 
traction rate t on S„, and velocity v on S„, the finite element 
algorithm presented in the previous chapter enables us to compute the 
(instantaneous) velocity v(x,t) and stress rate t(x,t). To h e more 
explicit, the information above suffices to compute the matrices [H], 
[G], {F}, and {P} on each element,* and hence, the velocity and stress 
rate throughout the body (see 7-^6 and 7-^7, and development thereof). 
The matrices [H] and {P} each depend upon the constitutive matrix W 
(through D) , and hence upon the quadrature point values of the stress, 
but in an inexplicit way.** 

We formally indicate the dependence of the nodal values of the 

velocity {v} » {v 1 ,v^, . . . , v N ^} (ND being the number of nodes) and the 

• _ • 1 *2 •G-, / 

quadrature point values of the stress rate (t) ■ tt ,t ,...,t } (G 


* * p. 

and if need be, the matrices [C], [C], [LJ, and (L }. 

** 

when the constitutive matrix ^ is a constant, the dependence 
of W on x is approximately affine; see Appendix B. 
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being the number of quadrature points) upon the nodal positions 
{x} ■ tx' ,x 2 , . . . ,x ND ), the quadrature point values of the stress 
{t } * {t ,t and the time dependent prescribed loads by 

•» ^ « -V 

writing* 

{v} * f[(x) »{ t) f tJ (3. I ) 


{t} - g[{x) , ( t } , t ] . 


( 8 . 2 ) 


Since each element node is associated with the same material point 
throughout a deformation, and likewise for each quadrature point, we 
may write each component of (x), (v), {t}, and { t } as 

x’ - X T (x',t) ; 

y 1 - X T (x',t) ; 

T 1 - 5 

t 1 - 0 /j!)f' • t (x 1 ,t) . 

T l ***' l 


The quadrature point values of the deformation gradient 


The functions f and g are introduced specifically as a 'short- 
hand' for the solutions of the finite element equations, as given by 
(7.1*6) and (7.1*7). from (7.1*6) and (7.1*7) it is clear that Integrations 
may be carried out on one element at a time. 
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- (v x T (x, t )) T 


X =x 


and its determinant 


J T * det 


depend on the nodal positions through the isoparametric shape functions. 
The reference configuration indicated by the subscript T is absolutely 
arbitrary and may be changed as frequently as desired. Typically it 
would be chosen on the basis of economy (e.g., to minimize storage 
requi rements) . 

The equations (8.1) and (8.2) may now be written as 


(x T > = f T [{x T Mt T },t] 


(8.3) 


(t T > * g T [(x T Mt T },t] 


(8.4) 


wh 


ere f and are defined by 


- fK{y- F T *t T ),'] 

T 


(8.5) 


g T [*,{t_},*]» [J T F“ 1 ]g[ F t * 


( 8 . 6 ) 


Equations (8.3) and (8.4) represent a system of nonlinear ordinary 
differential equations. On account of the complicated way in which 
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the finite element equations depend upon the stress, the functions f 
and (or even f and g) cannot be written explicitly. This fact makes 
stability analyses such as those of Cormeau [46] and Hughes and Taylor 
[20] Impossible for the equations (8.3) and (8.4). 

An Initial value problem may be set If initial values of X T 
and t^, the constitutive equation, and a program of loads are given. 

It is assumed that solutions of the initial value problem exist for 
sufficiently smooth and physically tenable Initial data, at least for 
some range of deformation from the Initial configuration. 

Numerical Integration of the Initial 
Value Problem 

The initial value problem posed by (8.3), (8.4) and appropriate 
Initial data is dependent upon the finite element equations as dis- 
cussed In the first section. From that same discussion, and from the 
presentation of the finite element equations (Chapter VII), it is also 
clear that the finite element- ini t ial value problem is predisposed to 
numerical integration. In this section we Indicate the types of 
numerical integration schemes suitable for the present problem, and 
mention a few important differences between the various types. 

The finite element- in it ial value problem may be integrated by 
single step explicit schemes, multistep explicit schemes, or (generally 
multistep) predictor-corrector schemes. A number of these schemes are 
discussed in the textbook of Conte and de Boor [47]. Three Important 
facts to be kept In mind when choosing a particular scheme are 

(1) the solution vector < CX T ( t N ) > , { t T ( t N ) > > at the time t * t^ 
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is of scalar dimension ND0F + 9*G, where NDOF is the number 
of kinematic degrees of freedom of the mesh and G is the 
total number of quadrature points. Storage required for 
implementation of different integration schemes can vary 
appreciably. 

(2) evaluation of <f ,g T > is expensive on account of the com- 
plexity of the finite element equations. 

(3) the functions f^ and g^ are generally discontinuous at points 
< (x T )» (t) > wh ich correspond to material yield surfaces. 

The multistep methods (explicit and pred i ctor-corrector) require 
relatively few evaluations of <f T ,g T > per step; this is an attractive 
feature. However, multistep methods are not self starting, the time 
step is not easily changed, they have relatively large storage require- 
ments (since several past values of <f T ,g T > must be carried along), 
and moreover, they cannot be expected to be accurate when the solution 
crosses a yield surface (since they are based on smooth polynomial 
interpolation of the solution over several time steps). 

On the other hand, the single step methods (explicit and 
predictor-corrector) are easily started, the time step size Is easily 
adjusted, and they have relatively small storage requirements. They 
can be expected to perform more favorably than the multi step methods 
when the solution crosses a yield surface since smoothing over several 
time steps Is not built in. The disadvantage of the single step 
methods is that a relatively larger number of evaluations of <f T »g T > 
are required per step to achieve a given accuracy when a yield surface 

is not crossed. However, the advantages of single step methods seem 
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to far outweigh the disadvantages. 

In the examples accompanying this work the Euler and classical 
second and fourth order Runge-Kutta methods were used. Details of 
these methods may be found in the textbook of Conte and de Boor [1*7]. 

We note that the second order Runge Kutta method Is equivalent to an 
Euler predictor and a single application of the trapezoid rule as a 
corrector. Errors of the Euler method were gauged (qualitatively) by 
step-halving and by comparison to results of second order integration 
for randomly selected time steps. Errors of the second order Runge- 
Kutta method were gauged in like manner; by step halving and comparison 
to results obtained by fourth order integration. The integration 
schemes used in the examples accompanying this work varied from problem 
to problem, and sometimes within a problem. Full details are given in 
the description of the individual problems, in the chapter following 
this one. 

We assert that the stress t_ (and hence t) integrated 
numerically satisfies LMB. As an example, consider the Eul er-Trapezoi d 
predictor-corrector pair for the stress at time (t N + h): 
predictor (Euler rule): 


{ ^° )(t N + h)) ■ { *r (t N )} +hl XT (t N )} 
{ ‘T° >(t N + h)l ' { iT (t N n + h( i T (, N )} 


(8.7) 


corrector (trapezoid rule): 
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tx' K+l) (v h )) - (x T (t N )) * ih(x T (t N ) +x{ K) (t N + h)> 

tt' K+,) (t N + h)} - (t T (t N )} + ih{t x (t H ) + t' K) (t N *h)) (8.8) 

To verify the satisfaction of LHB , we check to see that 

I 

/ t (K+1) : V„xdV - 0 (8.9) 

V ~ 

for all x(X) which vanishes on S (body force has been presumed to 

vanish, for simplicity). Elimination of t| K+ ^(t^ + h) from (3.9) by 

(K) 

use of the corrector (8.8), and assuming that (t N ) is balanced, 
yields 

JT [t T (t N ) + t^ K) (t N + h)| : V x xdV * 0 (?) (8.10) 


Since the stress rate is of the form 


t - V x $(t) , 


( 8 . 11 ) 


• • (k) 

the stress rates t, and t are of the forms 

~T ~T 


t ■ V x * ( t ) • t ^ « V 

lx X W ’ lx 


x $W(t + h) • 
X ~x VC N ' ’ 


( 8 . 12 ) 


(K) 

where $ and are defined by 

-T ~T ' 
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4tr 


X;i 


MV ’ (V X^T (t N )) *! ( V 5 

^ K> (t N + h) - (V x 4 K)(t N + h)) *! (t N + h) ' 

Elimination of t T and t^ from (8.10), and integration by parts 
affirms the satisfaction of (8.10). It is worthy of special note that 
the complementary work and energy principles provide no means whatever 
for checking the satisfaction of linear momentum balance by the stress, 
so it is of crucial importance that the numerical integration of the 
stress not introduce errors which lead to an unbalanced stress. This 
maintenance' of balanced stress, necessary In stress-based finite element 
algorithms, is the counterpart of maintenance of compatible deforma- 
tion, necessary in veioci ty-based finite element algorithms. 

The true stress is found, after t_(t M + h) and + are 

integrated, by the formula 

T(t N + hJ - (l/J T )(7 x x T (t N +h)) T • t x (t N + h) . ( 8 . 13 ) 

The true stress given by (8.13) necessarily satisfies LMB since 
t (t M + h) does. it does not appear to be possible to integrate the 
true stress explicitly with out causing it to become unbalanced, so 
no further consideration is given to that alternative. 

Angular momentum balance is satisfied only approximately by 
stresses computed in the present method. To keep the accumulated 
error small we embed the angular momentum balance as the stable 
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solution fl (t) - 0 of the following initial value problem: 


0 T + X« t - 0 0(0) - 0 ; 


(8.14) 


where the tensor fl is defined by 

—T 


9t 





(8.15) 


and X>0. In the course of numerical integration we adjust X at the 


beginning of the time step from t^ to t^ + ^ as 


X - 1/h = l/(t N+1 - t N ) . (8.16) 

The optimality of this choice is seen if (8. MO is replaced by the 
difference equation 


?t (t N + h) * (1 ’ Xh, 3?T ( V ' (8.17) 

Equation (8.14) accounts for the 'AMB residual' term in the finite 
element equation (7.24). 

Finally we note that frame indifference can be satisfied only 
In an approximate sense when one integrates the stress numerically. 

As an illustration consider the Euler method as applied by observers 
in frames which spin relative to each other. The first observer 
obtains 
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^T (t N +h)J " Ht (t N ) + h ! T (t N ) 


( 8 . 18 ) 


and the second observer obtains (for the same material point and the 
same time step) 

[t^(t N + h)J - t^(t N ) + ht|(t N ) . (8.19) 

We assume that t T (t N ), t T ( t N ) , ( t^) , t^(t N ), and Q(t) (the rotation 

between the two frames) are known exactly. According to the transfor- 
mation rules, the exact t^ and t^ satisfy, at each moment of time, 

t ' * t • Q T ; t' - t T * Q T + t_ • Q T (8.20) 

for arbitrary time dependent orthogonal Q. But at time (t N + h) 
equations ( 8 . 1 8 ) and (8.19) yield 

[t^(t N + h)] - [t T (t N + h)] • Q T (t N + h) - (8.21) 

[t T (t N )J • I($ T (t N ) +h§ T (t N )) - Q T (t N + h)] 

+ ht T (t N ) • [Q T (t N ) - Q T (t N +h)] . 

Since the right hand side of (8.21) does not vanish for all admissible 
Q and Q, the stress integrated by Euler's method will depend upon the 
frame of the observer. 
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Two courses of action are available. We might attempt to 
reformulate the tnitial value problem so that a frame Indifferent 
stress is Integrated, or we could attempt to integrate the stress in 
some special manner so as to remove the frame dependence. As pre- 
viously discussed, Integration of another stress besides t_ leaves us 
with no way to ascertain the satisfaction of LMB, so we disregard the 
first option above. 

For insight to the second option, let us consider integration 
of the stress when the deformation is homogeneous. Suppose that the 
spins o> and w' vanish in the frames of the two observers of ihe 
previous example. According to the transformation rule for spin 

i T • T 

u) = Q»a)»Q+Q*Q 

and since both and w 1 vanish, 

Q, » 0 ; Q ■ " constant . 

In words, all the frames in which the spin vanishes rotate as one with 
the principal axes of the deformation. Keeping In mind that the 
deformation is homogeneous, suppose we permit Euler's method to be 
applied only in frames in which the spin vanishes. Then (8.21) Is 


reduced to 



From this equation we surmise that a frame indifferent stress may be 
found In any frame If the classical Euler's method Is used only in 
those frames in which the spin vanishes. 

Let t T (t N + h) be the stress Integrated In the frame of the 
principal axes of the stretching (of the homogeneous deformation), and 
let + and t"(t^ + h) be the stress In two a r b i trary frames, 

determined from t (t..+h) as 

‘T' t N* h) 

Then the stresses t^(t N + h) and t”(t N + h) are related as 

t^(t N + h) - t^(t N + h) • Q" T (t N + h) • Q' (t N + h) . (8.24) 

Thus, the stresses a °d t^!(t N + h) are ^ rame independent. 

Equation (8.23) gives a clue as to how the stress may be 
Integrated in a general frame (when the deformation Is homogeneous), 
in (8.23) t (t.+h), the stress integrated In the frame of the 

~T N 

principal axes of stretching, is given by 

- t T (t H ) +ht T (t M > . 

which may be used along with the formulas (8.20) to write the first of 


*T ( V 


h) • Q 1 


(t N + h) 


(8.23) 


'x (t N + 


h) • O' 


(t N +h) 
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(8.23) In the form 



The prime is dropped in the equations below, but the discussion is for 
general frames. 

In a paper of related interest Rubinstein and Atluri [48] dis- 
cuss approximate solutions of the initial value problem posed by 
(8.26) for orthogonal Q.( t ) . In practice usually only w(t N ) is known. 
Then the best approximation for *Q(t^ + h) is given by 

Q T (t N ) • Q(t N + h) * | ^ sin(flh)a) N (8.27) 

+ 4r (1 - cos ( ft h) )w 2 
ft 2 ~ N 

*since only Q 1 T (t^) • q' ( t N + h) appears in (8.25), the initial 
value of q' only needs to be orthogonal, hence, Q. (t N ) ■ 1. 
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where 


ft 2 * i o> N : w N . ( 8 . 28 ) 

Since only an approximation for Q T (t N ) • Q.( t N + h ) is available, applica- 
tion of ( 8 . 25 ) through ( 8 . 28 ) as an integration scheme has the effect 
of integrating the stress in a frame in which the spin nearly vanishes. 
However LHB is precisely satisfied since the approximation for 
Q T (t N ) • Q( t N +h) is precisely orthogonal. 

Tor use in the finite element method one would replace by 
the mean spin (on each element) 

u - ~ / wdV . (8.29) 

N *V N 

Then we call the equations ( 8 . 25 ), ( 8 . 27 ), and ( 8 . 29 ) the 'modified 
Euler method'; they are summarized below: 

t T (t N + h) - [t T (t N ) + h(t T (t N ) + £ T ^ N )*w)] • Q T (t N ) • Q(t N + h) 

Q T (t w ) • Q(t w +h) ■ I - — sin(nh)uj + -i (1 - C os( h ) )oj 2 

N " N ~ n " a 2 

Q 2 - i u : u . (8.30) 

As long as the spatial mesh is fine enough to render the spin nearly 
constant on each element, this scheme is equivalent to application of 

the classical Euler's method in a frame in which the spin nearly 
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vanishes. 


By similar arguments we can establish frame indifferent Euler 
schemes for the true stress and Ki rchhof f stress as 

T(t N + h) - Q T (t N + h) • $(t )«[T(t M ) +ht*(t N )] • Q T (t N ) • Q(t N + h) 
where t*«*t + t»u)-u3*t 

and 

a T (t N + h) - Q T (t N + h) • Q(t N ) • + h a*(t N )] * Q T (t N ) * Q(t N + h) 

• JL • 

where o' ® a + a • uj - u) • a 

~T ~T ~T ~ ~T 

Of course Q^ft^) • Q(t^ + h) is computed approximately according to 
equation (8.30). 

To Illustrate the advantage of using the modified Euler method, 
consider the integration of the stress in a body which spins rigidly; 
that is 


t (t) - t • Q(t) ; 0(0) ■ I ; 

t T (t) “ t Q • Q(t); 

where t is a constant. An initial value problem for the stress t (t) 
~o ~T 

(in the frame in which body appears to spin) may be set as 


t_ (t) ■ -t (t) • 03 


-T 


-T 


-O 


t T ( 0 ) " t 

•»T ~Q 


where 03 q * -Q - (t) • Q(t) * constant. Euler’s method applied to this 
initial value problem gives the result: 
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t ( Nh ) « t • (I - hw ) N . 

The accumulated error grows without bound. On the other hand, the 
modified Euler method gives the exact answer: 


t T (Nh ) » t Q • 0( Nh ) 


1 1 2 

since Q(t) » 1 - rr s T n ( Q t)u ♦ — s- (1 - cos( !lt))u) , and, for constant 
~ - w ~o qZ — o 

ui , Q(Nh) - [Q(h)] N . 

In the integration of the examples accompanying this work the 
term (hu) was always so small as to make the classical and modified 
Euler methods Indistinguishable. However, in the general case one 
must take special precautions in the integration of tensors to insure 
that a frame dependence is not induced by the integration scheme. 

Stability of Numerical Integration of 
the Initial Value Problem 

It is possible that the difference between two supposed 
numerical solutions of a given Initial value problem is much larger 
than would be expected to arise from discretization error alone. As 
an example, consider integration of the stress in a material of the 
type (4.20) by the Euler method. We suppose, for the sake of 
simplicity, that e(t) Is given and I(t) ■ -2y(jyT'), so that the 
difference between two solutions satisfies 

ha* • [V(t + At) - V(t) ] : e(t) - (3 uy)At' . (8.31) 

•w m mr * — ** 
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If the elastic matrix and stretching are such that, in the Euclidean 
norm. 


1 [V(t + At) - V(x)3 : e(t)J/ |Ax|h o (8.32) 

*V ^ M •%* 

as ||At || 0, then for sufficiently small ||At || , equation ( 8 . 31) may be 

replaced by 


Aa* ■ -(3 vy)At‘ . 


(8.33) 


Defining Ao as Aa» J V? At' : At' , we may reduce (8.33) to a scalar 

O c +* ~~ 

equation in the invariant Ao: 


4r (Aa) - -(3uy)Aa . 


(8.3M 


For an initial value Ac(0) (small), the closed form solution of (8.3M 
i s 


Aa(t) 


Aa(0)e 


-(3UY)t 


(8.35) 


I 

Euler's method yields 


Aa N - Aa(0) (1 - 3UYh) N . 

It is clear from (8.35) that Aa decays to zero as time passes. 


(8.36) 


This 
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means that the closed form solution of the equation 

a* - V : e ( t ) + l ; t(0) - t (8.37) 

■** «s» % •»- O 

is stable with respect to sufficiently small perturbations of the 
deviatorlc part of T. On the other hand, the numerical solution (8.36) 
attenuates as time passes only if 

|(1 - 3PYh ) | < 1 . (8.38) 

This means that the numerical solution of (8.37) is stable with respect 
to small perturbations of the deviatorlc part of t' only so long as 
the time step h is bounded above as 


h | < 2/(3VY) 


(8.39) 


This bound is identical to the bound given by Cormeau [A6] (see 
equations 16 and 5** In this reference). It Is not surprising that 
time steps such as (8.39) are found to be necessary for stability of 
numerical solutions of the finite element-initial value problem pre- 
sented in this work. According to Hughes and Taylor [20]: 

The time step restriction of the Zienkiewicz- Cormeau algorithm 
... Is a stringent one in practice. For slowly varying loads, 
or when equilibrium response Is of prime interest, stability 
requires that time steps be selected which are much smaller than 
those necessary for accuracy. 

Argyris et al. [19] remark that this time step restriction amounts to 
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limiting the inelastic strain increment to be smaller than the elastic 
strain. Since the elastic strain is usually very small in metals such 
as those used In structures, this implies that a finite deformation 
analysis would entail an Impract i ca 1 ly large number of steps. 

The work of Kanchi et al. [49] and Atluri and Murakawa [14] 
suggests the modification we now describe. To improve the estimate of 

p 

the inelastic strain increment in a time step, we replace e (t^) by an 
estimate of the mean value of the inelastic stretching in that time 
step: 


J P( I ( V 


eh)) i e P (T(t N )) +eh 


a* 




(8.40) 


where the parameter 0, 0<8^1, serves to locate the time at which the 
mean value is achieved. As 0 goes from zero to one, the estimate of 
the creep-stretch becomes increasingly more conservative. 

Equation ( 8 . 40 ) may be introduced to the finite element 
algorithm through the constitutive equation; (4.20) becomes 




i * Ee 


where 





+ 0h 




From V Q (8.41) we derive W a just as we derived W from V: 

2: s 


(8.41) 
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y 

as 


e 


>9 



J jkl 


*< T ik 5 lj 


5 ik T i j 5 


h m ‘He : £ P " ” { ie + P : £ P • (8.1*2) 

When a material which exhibits relaxation is to be analyzed, W- and Z Q 

are Introduced to the finite element algorithm for W and Z. 

* - 

If Euler's method is used to integrate the Initial value problem 
which results from the modification described above, then, in the 
terminology of Argyris et al. [193, an explicit 'forward gradient 
scheme' is recovered. If the 'gradient' (de /dt) is evaluated at 
T» Tg ■ (1 - 9)t n +0t n+ | , then an implicit forward gradient scheme, 
counterpart to that proposed by Hughes and Taylor [20] is recovered. 

p 

Finally, if the 'gradient' (de /dT) is replaced altogether by a 

■M 

function g ft , defined through 


£ P (t) - g(t) : T 


fe “ I { Ie ) ; 


(8.1*3) 


then an implicit ‘finite approximation technique,' counterpart to that 
proposed by Argyris et al. [193 Is recovered. The implicit schemes 
must be solved by Iterating on each time step, keeping 0 fixed. The 
iterative schemes amount to special predictor-corrector techniques. 

An Important fact that is exploited in our numerical studies Is that 
for the material ( 8 . 31 ), with yy constant, all of these schemes are 
equivalent. 
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CHAPTER IX 


EXAMPLES: FINITE DEFORMATION PROBLEMS 

Int roduct I on 

In this chapter we present several examples as demonstrations 
of the feasibility and performance of the finite element algorithm. 

All of the examples may be described either as plane or axi symmet r i c, 
so we begin by discussing specializations of the algorithm for such 
prob 1 ems . 

The examples fall into two categor ies--homogeneous deformations 
and inhomogeneous deformations. By treatment of homogeneous deforma- 
tions (for which closed form solutions are known), several important 
aspects of the performance of the finite element algorithm can be 
clearly identified and studied. The studies of inhomogeneous deforma- 
tions, the results of which are compared to both analytical and other 
numerical results, indicate the potential of the algorithm for treat- 
ment of problems of technological interest. 

In the discussions of the examples the time integration scheme 
used is indicated as Euler, Runge-Kutta second order (RK2), or Runge- 
Kutta fourtn oraer (RKm). In all cases the classical schemes (i.e., 
those described in [1*7]) have been used. 

Plane Strain 

All but one of the deformations studied in this chapter are 
plane strain in character. Just as for formulations using ordinary 
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stresses, a number of the components of the velocity, spin, and stress 

rate vanish if a Cartesian coordinate system is chosen with one axis 

2 

normal to the plane of deformation. We have chosen the x coordinate 
line to be normal to the plane of deformation, so that the velocity, 
spin, stress rate, and stress are of the forms 


v ■ v^e^ + v^e^ 


13 ^ 31 

co *» a) + ur e^e^ 


• «1 1 • 1 3 

t » t e^ + t e j 


x ; 22 
+ 1 - 2-2 


•31 *13 

+ t e^e 1 + t e^e^ 


T * T 1 1 e 1 e 1 + T 1 ^e^ 


+T s 2 s 2 


31 33 

+ T - 3-1 + T - 3-3 


None of the components depends upon x . The velocity is represented 
on each element as 


NQ 


Es. 

i-1 
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The shape functions N. are described In Appendix C. Similarly the spin 
and stress rate are represented as 

NW NT 

“ ■ E ■ i - E aji 6 ' 

1*1 i»l 

and those shape functions are given In Appendix C also. 

We note that this approach requires minimal specialization in 
programming for the particular case of plane strain. The plane strain 
condition is not satisfied a priori; that is, 

Se 22 * ^-2-2^ : w 1 : * 0 

for arbitrary 6r. Rather, E 22 * ^ ^°^ 1ows f rom the stationary con- 
dition (a component of 7.1): 

f [-c 2 l (t. S )] «J 22 dv - 0 . 

In using the finite element algorithm the plane strain condition is 

only satisfied approximately. In practice a qualitative check for 

satisfaction of the plane strain condition can be made by seeing that 

22 

the stress component T and the mean stress are nearly equal. This 
method for checking * 0 works so long as the inelastic stretching 
is proportional to the stress deviator (in the constitutive equation). 

As an alternative to approximate satisfaction of the plane 
strain condition, it is possible to 'split' the constitutive equation 
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Into two equations. The first Involves only the In-plane components of 

the stretching and stress rate, all that Is actually required for the 

finite element algorithm for plane strain problems. After In-plane 

components of the stress rate and spin have been found, the stress rate 
•22 

component t can be assigned so as to give precise satisfaction of the 
plane strain condition. This alternative is attractive from the point 
of view of efficiency, for storage could be reduced and the analytic 
Inverse for the constitutive equation could be used, but it requires 
considerable specialization in programming. Because of the inflex- 
ibility of this approach, it was not pursued in programming. 

Axi symmetric Deformations 

One of the deformations studied In this chapter is axisymmetric 

in character, that of creep of a pipe from Internal pressure. We have 

used a right circular cylindrical coordinate system to describe the 

1 2 3 

problem, indexing the coordinates as x * r, x * 0 , and x *z. Of 
course the z axis is along the centerline of the pipe, and r is a 
constant on the interior and exterior surfaces of the pipe at any 
particular time. The velocity, spin, stress rate, and stress are of 
the same forms as for plane strain; that Is 

v ■ v^ej + v ^?3 

13 A 31 

u - uj 5,53 + <*> S3S1 

• *11 *13 *22 *31 *33 

t - t e j e j + t ^§,53 ♦ t S2S2 + * -3-1 + * -3-3 
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T " s l-1 


13 22 31 33 

+ T + T e 2 e 2 + r 5 e^e 1 + 


2 

None of the components above depends upon x ■ 0. The representations 
for the velocity, spin, and stress rate are of the same forms as for 
plane strain; the shape functions are described in Appendix C. 


Homogeneous Deformations 

Through the study of homogeneous deformations various important 
•aspects of the performance of the finite element algorithm can be 
identified and studied. Since closed form solutions to problems of 
homogeneous deformation are widely available, questions of the accuracy 
of the finite element solutions can be resolved quickly and absolutely. 
If we immediately engaged problems complicated by inhomogeneous 
deformation, the accuracy of any solution we obtained would be no 
more than a subject for speculation. It was demonstrated by the 
example In Chapter VIII that homogeneous deformations are as difficult 
to integrate from the point of view of time step stability as 
Inhomogeneous deformations (since the same time step bound was found). 
Thus, homogeneous deformations are convenient subjects for studies of 
time step stability, and in the present case, for studies of the 
effect of the forward gradient scheme on accuracy. Finally, the results 
of this study serve to underscore the fact that the material models 
themselves are too idealized to be used in problems of technological 
Interest when strains are very large. 

Finite Plane Extension 

We begin our study of homogeneous deformations by considering 
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finite plane extension of (i) the hypoelastlc material (4.30, 

(ii) a hypoe last! c/p last Ic material, and (111) a creeping viscoplastic 
material. The geometry of the specimen for these examples is given 
In Figure T. These examples serve to demonstrate the relative 
efficiencies (accuracy/cost) of the Euler, RK2, and RK4 time integra- 
tion schemes, the performance of those schemes when the loading path 
crosses a material yield surface, and the effect on accuracy of the 
forward gradient scheme. 

Hypoelastic Material . For the first problem we consider the 
hypoelastic material defined by the constitutive equation (4.31): 

a* » A(l : e) I + 2 ye . (9-0 

«V V •* ^ 

It is convenient to introduce the normalized stress s, defined as 

"V 

s*x/2y. Then (9-0 may be written 

!* - (ttj*) |i 1 S’! ♦ « • <9 - 2) 

For homogeneous plane extension of this material from a stress free 
state, the stresses are given by 

s n (l) - 0 

s 22 (!) - vs 33 (l) (9.3) 
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where v» (1 - 2v)/(1 - v) . In Figures 2, 3, and 4 the stresses found 
by application of the finite element algorithm are compared to the 
closed form solution (9.3) for v-1/3. Euler, RK2, and RK4 time stepping 
algorithms are used. The Euler algorithm underestimates the strain- 
softening of the material slightly, but the RK2 and RK4 algorithms pro- 
vide data indistinguishable from the closed form solution. The reader 
should note that the stretch increment for each time step is (0.01), 
(0.02), and (0.04), for the Euler, RK2, and RK4 algorithms, respectively. 
Thus, the computational effort is the same in each of these three 
cases (1 element stiffness matrix evaluation per stretch Increment of 
0.01). In view of the differences In accuracy, we rank the RK4 
algorithm as most efficient, followed by RK2 and Euler. 

One of the methods used to get qualitative estimates of the 
local error of the Euler and RK2 algorithms is based on the assumption 
that greater accuracy will always be achieved by the next higher order 
time stepping algorithm when the solution is smooth. To check the 
error of the Euler algorithm on a given time step, we integrate the 
time step a second time using the RK2 algorithm. The local error is 
then assumed to be of the same order of magnitude as the difference 
between those two solutions. The error of the RK2 algorithm may be 
checked in a similar manner using the RK4 algorithm as a 'reference.' 

Plastic Material ♦ As a second example we consider plane 
extension of a hypoelastic/plastic material from a stress free state. 

The specimen geometry is Identical to that of the previous example 
(see Figure 1). The material is characterized by uniaxial test data 
as 
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Figure 2. Stress 
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Figure 3. St 
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Figure *». Stress Accompanying Plane Extension of a Hypoelastic Material — RK 1 ! Integration 





(9.M 


• l 1 

a o /e 1 1 * J o E 


(T 11 < T y ) 


J o' /e U ’ J o E 




1 UN-1 ) -1 

y 


, 11 ^ v 

(T ^ T y ) 


The associated constitutive equation is 


e ■ e e + e P - V : a* 


(9.5) 


) 0* - ( 

t) (' :°*)l ) 

J - \ 

E/ * * ' ( 


. I 

(\ 1') 

(T : 0*) T ) 

\ 2 / 
I*. 44) h 

is defined as 

- 1/E * 

(1/E) | N 


N-l 


- 1 


(9.6) 


(9.7) 


where I ' * t' :t', and a* 0 or a* 1 as the material is elastic or 
plastic. It is convenient in this case to normalize the stress by its 
yield value in plane extension: 


s - ST t/2t 


(9.8) 


The stress accompanying plane extension from a stress free state of 
the material (9.5) are easily found when v*Jas 
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(9.9) 


s 1 1 (1 ) - 0 ; s 22 (1) = i s 33 (l) 

s 33 (l) - [ln(1)/ln(!)J s < 1 

s 33 (1) = [ln(l)/ln(l y )] N s> 1 

where 1 is the nominal strain at the initial yield of the incompres- 
sible material in plane extension: 

In ( l y ) * 7T T y /2E . (9-10) 

In the examples we take (E/t ) * 200 , N»4, andv«*(l/3). In Figures 
5, £>, and 7 the stresses found by application of the finite element 
algorithm are compared to the stresses in the incompressible material 
with the same Young's modulus (E/t y ) and hardening exponent N. The 
discrepancy in the elastic range is due entirely to the difference in 
the Poisson ratio; for both v * (£) and v = (y) the stress T *= vr , 
as we have determined it should in the previous example (see 9*3). We 
note that the numerical solution 'overshoots' the yield surface for 
Euler, RK2, and KKk algorithms, but the error appears to be least for 
the higher order methods. As plastic deformation progresses, the 
stress in the compressible body falls slightly below the stress in the 
incompressible body. This is to be expected, for the compressible 
body is more compliant. One other point made by this example is that 
the deformation in the plastic range (beyond about 2% stretch) is 
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Figure 5. 


Stress Acccwnpaoy i ng Plane Extension of a Plastic Material — Euler Integration 



STRETCH 


Figure 6. Stress Accompanying Plane Extension of a Plastic Material — RK2 Integration 
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Figure 7. Stress Accompanying Plane Extension 


(IP I' 'f M 




essentially isochoric. This is not because the material is any less 
compressible in the plastic range; rather, it occurs because of the 
drastic lowering of the shear rigidity of the material by the plasticity 
mechanism. The condition number of the constitutive matrix V can be 
used as an indicator of the relative shear rigidity; as that number 
decreases, an unconstrained deformation becomes more and more like the 
deformation of an incompressible body. 

In Figure 8 of the stress accompanying the deformation out to a 
stretch of 1.92 is plotted for the same materials as above. The RKA 
algorithm was used for this integration, with a stretch increment of 
(0.04) . 

Creeping Viscoplastic Material . As final examples of plane 
extension we consider a viscoplastic material which creeps (i) from an 
initially stressed state, and (ii) from a stress free state. The 
specimen geometry is the same as that of the previous two examples. 

The material behavior is governed by 


e 




y : 9l) 


(9.11) 


where s»t/2ij, s*=o :V /2y, and T=l/3uy. When the material is incom- 
pressible, the stress accompanying plane extension are easily found as 


s 22 (t) 


is 33 (t) 
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Figure 8. Stress Accompanying Plane Extension of a Plastic Material — to 1 = 1.92 


I j 


I 


ill 


s 33 ( t ) - s 33 (0)e' t/T 


(9.12) 


+ 2 r * £ 33 (0e -( t "0 /T dc . 

0 

For the stretch history 1(t)*l+Vt, the stretching e 33 (t) is 

e 33 (t) • V/(l + Vt) . (9.13) 

The stresses for this stretch history take the form 
s n d)-0 
S 22 (l) - *s 33 (l) 

s 33 (1) - S 33 (0)e‘ (1 ' 1)/L + 2e‘ 1/L |ei (1/L) -Ei(1/L)J (9.U) 

where L**VT, and Ei(’) is the exponential integral, defined 
Ei(x) » J ( e z /z)dz 

-CO 

Values of Ei may be found in tables [53]* In the present case Ei was 
evaluated by a subroutine in the IMS L Mathematical and Statistical 
Library (FORTRAN). The subroutine (MMDEI) returns a value to the call- 
ing program which is based on interpolation of tables. 

-1 A -1 

In both of the cases which follow we take V»10 sec 

(nominal stretch rate) and T»* (^) x sec. Cormeau's [A6] time 

step bound is h < h « 2T. 

& 
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In the first case we assign the Initial stress as if it arose 


from viscosity alone; 

s n (0) - 0 ; s 22 (0) - *s 33 (0) - VT . 


We take time steps h ■ (90/56 ) h , 2 (90/56 ) h , and ^ (90/56 ) h (for 

c c c 

stretch increments of 0,01, 0.02, and Q.0*», respectively), applying 
Idler, RK2, and RK*» algorithms, respectively. The stability parameter 
9 is set as 9»i. The stresses found by application of the finite 
element algorithm are compared to the closed form solution (9*1*0 for 
the Incompressible body with the same characteristic time T in Figures 
9, 10, and 11. It is apparent from these figures that the numerically 
integrated stress is slightly greater than the stress in the incom- 
pressible material, which is surprising. One would expect the more 
compliant material to have the lower stress. 

We note that the introduction of the stability parameter 9 has 
the net effect of replacing (9.1 1) by 



Dividing through this equation by , setting V"i, and 

defining T„ * T + '•'h oive? 

v 3 
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Figure 9. Stress Accompanying Plane Extension of a Viscoplastic Material--Euler Integration 
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Figure 10. Stress Accompanying Plane Extension of a Vi$coplastic Material--RK2 Integration 
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(r) s - i* - {rrh) U : ;*>i * r(: - y ( ; '• i>i) • «-'5) 


But this equation has already been Integrated for plane extension. 
From (9.12) we can Immediately write 


* n (t) - 0 ; s 22 (t) - is 33 (t) 


.»(t> - ♦ f/' «»(c).* (, ‘ 5,/T# 


dc . (9-16) 


Now we consider the special case that £ J,? (t) is a constant. Integra- 
tion of (9.12) gives 


s 33 (t) - s 33 (0)e” t/T + 2Te 33 (0) (1 - e' t/T ) . (9.17) 


At times t much later than t * T the stress attains a steady state 
value of 2e 33 (0)T. Integration of (9.16) yields 

.»(«) - s”(0).' t/Te ♦ lft»(0)(l - «‘ t/Te ). (9.18) 

and at late times Sg 3 (t) 2 e 33 (0)T. Thus, the steady state value 
of the stress in plane extension is unaffected by the introduction of 
the stability parameter 8. The important difference between (9.17) 
and (9*l8) is the rate at which the stress approaches the steady state 
value. It is clear that s 33 (t) and ( t ) are related as 
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33 (,(,♦»)) . 

Simply put, the stress Sg^ 'lags' the correct stress. This Is exactly 
what Is. observed In Figures 9, 10, and 11. That lag becomes Increas- 
ingly apparent as the order of the Integration scheme is increased, 
since the numerical solution tends towards Sg^(t), not s^(t). We 
conclude the discussion of this example by remarking that this lagging 
of the stress must be expected any time the stability parameter 0/0, 
and the lag increases with (8h/T). 

In the second case of plane stretching of a creeping viscoplas- 
tic material we set the initial stresses to zero. As before, 

V*10 ^sec ' (nominal stretch rate), T * i { 56/90 ) 10 1 2 sec (character i st i c 
time), and Cormeau's [**6] time step bound is h<h * 2T. The stress in 
this example differs from the s ess in the previous example solely 
because of the different initial stress, and that difference 
attenuates like e . For times t much later than t * T, the stress 
In this example is indistinguishable from that of the previous example. 
We use Euler's method for integration, taking time steps of 

h»i(90/56)h , (90/56 ) h , and 2(90/56)h (for stretch increments of 
c c c 

0.005, 0.01, and 0.02, respectively). We set the stability parameter 
e to unity, ?=? ( g i •. i ' 5 (6H/T- ( 90 / 56 ) , 2(90/56), and 4(90/56), 
respectively). In Figures 12, 13, and 14 the stresses found by 
application of the finite element algorithm are compared to the 
closed form solution (9-12) for the incompressible body with the same 
characteristic time T, and also to the closed form solution for the 
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Figure 13* Stress Accompanying Plane Extension of a Viscoplastic Materi a 1-- (0h/T) = 2(90/56) 
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Figure IV. Stress Accompanying Plane Extension of a Viscoplastic Material-- (Oh/T) = Vi (90/56) 




It is 


incompressible body with the 9-modiFied characteristic time Tg. 
clear from these figures that only a qualitative estimate of the 
transient stress response is given by the finite element algorithm; 
and that estimate is degraded as (0h/T) Increases. More Important 
though is the fact that use of a higher order time stepping algorithm 
(for the same (0h/T)) cannot improve the accuracy of the numerical 
solution, since that solution would be drawn closer to Sg, not s. It 
is therefore senseless to use methods other than Euler's when 9h is of 
the same order of magnitude as the characteristic time T. 

Finite Rectilinear Shear 

We continue our study of homogeneous deformations by considering 
finite rectilinear shear of (i) the hypoelastic material (9.2), (ii) a 
second hypoelastic material which resembles an e last ic-perfect ly plas- 
tic material, and (iii) a creeping viscoplastic material. The geometry 
of the specimen is given in Figure 15. These examples serve not only 
to further demonstrate the performance of the finite element algorithm, 
but also to portray aspects of the finite deformation behavior of the 
materials themselves. Consistent with the conclusions of the pre- 
vious section, we use only higher order integration schemes for the 
two hypoelastic materials (on the basis of efficiency). It happens 
that for this completely constrained deformation stable integration of 
the stress in the viscoplastic material may be achieved without 
relying on the forward gradient technique. We take advantage of this 
situation by comparing the accuracy of the finite element algorithm 
with and without the forward gradient scheme for the same deformation. 

From the results one can only conclude that further research is needed. 
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Figure 1£. Rectilinear Shear Specimen 
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Hypoelastic Material . We consider for a second time the 


hypoelastic material defined by (9,2). For homogeneous rectilinear 

1 2 

shearing of this material from a stress free state, v «v "0, 

3 • 1 

v ■ ex , the normalized stress s is given by 


s^(e) ■ -s 11 ^) - i(l - cos (e) ) 

•** lit- 

s'^(e) » i sin(e) 


(9.19) 


We compare the stresses found by application of the finite element 
algorithm to the closed form solution in Figure 16. The RK2 
algorithm was used, with time steps corresponding to (nominal) shear 
increments of (0.16). As discussed in Chapter IV, constitutive 
equation (9.2) is invalid beyond the nominal shear strain e*iir. 

This is of little consequence if one is interested only in metals such 
as used in structures, because some mechanism of inelasticity always 
sets in long before such a large shear strain as e ■ i it is reached. 
However, since (9*2) Is typically used to njpdel the elastic part of 
the stretching in constitutive equations for inelasticity, one must 
ask whether or not similar periodic behavior will be observed when a 
hypoelastic/plastic or hypoel ast i c/vi scopl ast i c material is subjected 
to such large shear strains. 

A Second Hypoelastic Material . We consider a second hypo- 
elastic material whose constitutive equation resembles that of an 
elastic-perfectly plastic material. This material is studied in 
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Figure 16. Stress Accompanying Rectilinear Shear of a Hypoelastlc Material 
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place of a truly plastic material because of the availability of a 
closed form solution [27]* It would have been possible to 'splice 1 
solutions for an elastic material and the present material at the 
yield surface to obtain a closed form solution for a truly elastic- 
plastic material, but there would have been little to gain in the way 
of illustration. The new hypoelastic material is defined by 

a* - 2 ye + X ( I : e) I - ~ - T - (r * : e)r' . (9.20) 

| t 2 ~ 

3 y 

By defining K « y (t^/2p) , introducing the normalized stress 
s a T/2p, and prescribing rectilinear shearing, we deduce 

s * « e - -V s(s : e) . (9.21) 

For materials such as those used in structures the parameter K is 
typically of the order of (0.01). The stresses accompanying rectilinear 
shearing from a stress free state were found in [27] as 


s 13 (4>) - i Ksin(<J>)/(i + K 2 )* 

s 3 3(<j>) - -s n (<j)) - iK 2 (l - cos (<p))/(i + K 2 ) 


( 9 . 22 ) 


where 


tan(i$) 



(±z 

2K 



K 2 < i 
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In Figure 1? we compare the stresses found by application of the finite 
element algorithm to the closed form solution (9 - 22 ) or the special 
value K* (0.10). The RK2 time stepping algorithm was used for the 
computation. The closed form solution has a peak at [27] 

e - - — r In i. * ( * ~ ^) — - 0. 656*131 

(i - K 2 )* 

but this is completely lost in the numerical solution. At shear strains 

larger than e above, the shear stress s^ is a decreasing function of 
o 

e; one must conclude that the model (9-20) is unacceptable beyond that 
level of strain. 

Creeping VTscoplastic Material . As final examples of recti- 
linear shear we consider a viscoplastic material wh creeps from a 
stress free state at several different rates. For a second time we 
use the material defined by (9.11). The specimen geometry is the same 
as in the previous two examples. The stress accompanying the deforma- 
tion Is 

s^ 3 (e) ■ -s 1 1 (« ) » s ] J\ E-exp (-~)[sin(e) + Ecos(e)]j 

s 1 ^ (e) » spjl-exp (-jjp-)[cos (e) - E s I n (e ) ] j (9.21) 

where E«eT, and s^ ■ iE/(l + E^) . The characteristic time T is the 
same in this case as in the plane extension examples. We took values 

of E as (0.5), O.0), and (2.0), corresponding to nominal shear strain 
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Figure 17- Stress Accompanying Rectilinear Shear of a Second Hypoelastic Material 




rates of (45/28) x l(f 12 , (1*5/1 A) x 10* l2 , and (1*5/7) x 10~ 12 . respectively. 
The reader should note that these strain rates are indeed very slow. 

The stresses found by application of the finite element algorithm are 
compared to the closed form solution in Figures 18, 19, and 20. The RK2 
algorithm was used for time integration, with time steps corresponding 
to nominal shear strain increments of (0.2). The first stress peak 
is at e*iTr, and the final shear strain is (8.0). The reader should 
note that the behavior is (qualitatively) similar to the purely 
hypoelastic behavior (period) and to rigid-viscoplastic (or viscous 
fluid) behavior at late times. 

In the three cases above it was not necessary to use the forward 
gradient scheme since the deformation was completely constrained. How- 
ever we can still introduce the stability parameter for the pu-oose of 
finding out what error it leads to. Just as for the examples of plane 
extension, it is possible to integrate the modified equations in 
closed form. The solution of the 0-modi fied problem is found by 

replacing E by E(l+0h/T) in (9.23) everywhere except in the numerator 
1 3 

of s . The stresses found by application of the finite element 
algorithm, with 0*1, are plotted in Figure 21, along with the 
solutions of the true and modified problems, for the case E*2. All 
the other data is unchanged. The reader should note that at late 
times the two solutions do not coincide; rather, we get the ratios 


Sg 3 /s 13 - (1 + E 2 )/(l + Eg) 

(9.2A) 

11/ 11 ^ ft /r w 13, 13> 
s Q /s -*• (Eg/E) (sq /s ) . 
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Figure 18. Stress Accompanying Rectilinea- "hear of a Viscoplastic Material — E = 0.5 
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Figure 20. Stress Accompanying Rectilinear Shear of a Viscoplastic Material — E = 2.0 
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From (9-15) it is easily seen that the influence of the parameter 6 
will vanish at late times (for any deformation) only so long as s* 
vanishes at late times. In other words, the forward gradient scheme 
presented here will only give correct steady state values of the stress 
for deformations in which the spin or the total stress vanishes at late 
times. Rectilinear shearing does not meet either of these conditions. 
Further research is needed to determine whether or not this defect in 
the forward gradient scheme, which is probably only important for large 
deformations, can be overcome. 

Inhomogeneous Deformations 

We conclude this chapter with three studies of inhomogeneous 
deformations. In the first, a study of creep In a pipe due r o internal 
pressure, the numerical solution is compared to closed form solutions 
for elastic and rigid-plastic bodies (at different stages of the 
deformation), as well as to other numerical results. In the second 
problem, a study of the onset of necking in the plane tensile test, 
results are compared to closed form results for an incompressible body 
as well as to other numerical results. In the third problem, a study 
of plane void growth in a creeping viscoplastic medium, comparison is 
made to other numerical results. As shall be seen, the stress-rate 
based finite element algorithm compares very well In every case. 

Creep of a Pipe from Internal Pressure 

The problem geometry and boundary conditions are given In 
Figure 22. The problem may be analyzed in three stages: (1) rapid 

elastic inflation, ( 2 ) transient stress stage, and ( 3 ) large strain 
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creep stage. The first two parts of this problem have been studied by 
Greenbaum and Rubenstein [5*0, and Plan and Lee [551- The creep 
behavior is characterized by the uniaxial relation 


e 


P 

11 


Y (T n )V'-' 


(9. 25)* 


where Y=* 2.073 x 10 ^ (hr) ' (N/mm^ ) and M * 0. 7. The elastic 

and inelastic components of the stretching are 


■•■( if )!* - ( t ) ; ?*>2 »•») 

~ P * r Y ( Vi i' : i' ) (Mt M ~ 1 )x' 

As in [551, time integration is carried out with respect to the 

M 1 /M 

parameter X*t . Real time is recovered after integration as t ■ A . 

In the example we take E * 1 . 379 xlO^Pa, \) = 0 . 45 , 9*1, and use Euler 

integration exclusively. 

Rapid Elastic Inflation . The pipe is taken from a stress free 
state to an elastically strained state in a single Euler step by 
imposing a nominal pressure rate on the inside wall of the pipe of 
T, - 2517 (iviiim • sec j lor a A-time step of 0.001 (sec) . The pres- 
sure in the pipe at the end of elastic inflation Is 2.517 (N/mm ). 

In Figure 23 the stresses found by appl ication of the finite element 


when the parameter M< 1, the material age-hardens . 
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Figure 23. Stress Profiles in Pipe Crfcep Specimen 
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algorithm are compared to those predicted by the infinitesimal 
displacement theory of elasticity. As is evident, the two solutions 
are virtually indistinguishable. 

Transient Stress Stage . In this stage the creep mechanism 
causes the stresses to redistribute themselves through the thickness 
of the pipe, and the pipe itself to expand slightly. In this stage 
the nominal traction rate on the inner wall was set to zero, and the 
A-time step was set as AX* 0.5 (sec) . In Figure 23 the new distribu- 
tion of stresses is compared to the stresses that would be obtained if 
the material were rigid plastic. These latter stress distributions 
are given by Hult [56]. Again, the two solutions are virtually 
indistinguishable. In Figure 2A the small deformation displacement 
history of the inner wall of the pipe Is plotted against the result of 
Pian and Lee [55]. 

Large Strain Creep Stage . This stage is a continuation of the 
relatively steady creep that characterized the latter part of the 
previous stage. The only difference In the computation was the A-time 
step was set to A00. In Figure 25 the maximum hoop stress history 
(at the outside wall) is plotted against the history for the rigid- 
plastic material. In this stage it Is Important to note that the 
traction boundary condition on the inside wall, T^ - 0 does not corres- 
pond to constant pressure; rather fnom the equation (3.36) we find that 
the pressure rate is 

p ■ - (a/a)p (9.27) 
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Figure 24. Small Displacement History of Inner Wall of Pipe 
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where a is the instantaneous inner radius of the pipe. Integration 
gives 


P - (a Q p o )/a . (9.28) 

This traction boundary condition was chosen simply for convenience. 

For constant pressure, the nominal tract i on" rate and velocity on the 
inner wall would have” bo satisfy 

Tj - (a/a)p Q - 0 . (9.29) 

This boundary condition could be dealt with by iteration, using (9.29) 
to form 'residual loads,' or by recasting the traction boundary con- 
dition (7-2) as 

/ $ [t, - (a/a)p]Sa d5 « 0 . (9-30) 

a J 

The Onset of Necking in Plane Extension 

ft is well known that in tensile tests of metals (plane or 
uniaxial) the deformation prior to the attainment of the maximum load 
is essentially Homogeneous. At some point after the maximum load 
point a 'neck' forms and the specimen fails almost immediately. In 
experiment, the point at which necking begins is sensitive to the 
slenderness of the specimen, the rate of strain hardening of the 
material, and the presence of geometric imperfections or inclusions. 

Mathematical ly It Is possible to consider specimens absolutely 
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free from geometric imperfections, and throughout which material 
properties are perfectly homogeneous. For such a ’perfect' specimen 
homogeneous stretching (such as discussed earlier In this chapter) to 
any extension is a solution of the general boundary value problem. 
However, to reach configurations of pure extension far beyond the 
maximum load point, the perfect specimen must pass through a (possibly 
infinite) number of configurations from which bifurcation is possible. 
For the classical elastic/plastic solid with a smooth yield surface 
the first possible mode of bifurcation is that of necking. 

Within the past few years mathematical analyses of bifurcation 
from configurations of pure extension have been presented by Hutchinson 
and Miles [57], Miles [58], and Hill and Hutchinson [59]. The first 
two of these are concerned with the onset of necking in cylinderical 
and rectangular specimens of an incompressible el ast i c-pjast i c material. 
In the third an extensive study of general bifurcation phenomena of 
incompressible materials in the plane tension test is presented. For 
the classical elastic/plastic solid this study indicates that the 
first possible mode of bifurcation is that of necking. 

Here we present a numerical study of the onset of necking of 
(hypo-) elastic/plastic solids in the plane tension test. To aid 
in the comparison of our results to those of Hill and Hutchinson [59], 
certain of their notations are adopted; these are explained as they 
are presented. The material we consider is slightly compressible*; 


* 

the loading surface does not depend upon the pressure. 
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in order to compare our results to those of Hill and Hutchinson we 
must assume that the effect of this compressibility is slight, though 
not necessarily ignorable. So that the performance of the present 
finite element method may be contrasted to the performance of a 
velocity-based finite element algorithm, we have chosen particular 
materials identical to those used in the bifurcation study of Burke 
and Nix [60]. This appears to be the only other numerical study of 
bifurcation of classical elastic/p last tc materials in plane extension 
in the literature. 

Finally, we investigate the sensitivity of our results to 
variations in the number, shape, and type* of elements in the finite 
element mesh. This latter study serves not only to graphically 
demonstrate the stability of the finite element algorithm, but also to 
help characterize the approximation thus obtained. 

Bifurcation Analysis . As discussed in Chapter VI, a condition 
sufficient for uniqueness of a deformation of an (hypo-) elastic/ 
plastic body is that the second complementary virtual work principle 
have only the trivial solution for the 'linear comparison solid' when 
homogeneous boundary conditions are imposed. In as much as the finite 
element algorithm is based on that work principle, we equate the 
uniqueness criterion with the condition that the finite element 
equations have only the trivial solution for the linear comparison 
solid when homogeneous boundary data is Imposed. Searching for possible 


the type is determined by the number of boundary nodes, the 
number of spin parameters, and the number of stress parameters. 
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bifurcation points amounts to looking for configurations of the body 

in which the global stiffness matrix [K*] is singular. The subscript 

L is to indicate that the stiffness matrix is formed as if loading 

occurs throughout the plastically stressed portions of the body; in 

extension problems this condition can be satisfied a posteriori by 

judicious superposition of the homogeneous and necking modes. 

We consider a specimen of initial length 2a^ and thickness 

2a^. We assume that the bifurcation mode will be symmetric in the 

sense that the velocity field may be reflected across the x^ axis 

(see Figure 26). This is consistent with usage of the adjective 

'symmetric' by Hill and Hutchinson. The finite element mesh is over 
1 3 

the area a x 2a . The specimen is composed of the hypoel as t i c/p last i c 
material of equations (9-5) through (9.7), and the subsequent dis- 
cussions pertain only to that material. 

To ease the comparison of our results to those of Hill and 
Hutchinson [59] we introduce the plane strain tangent modulus 
defined by 

. (9.31) 

o 33 

Assumi^- :h = ‘ s' V u 1 ' c jrcat ion point is approached the stresses 
satisfy (approximately--see Figure 8) 


22 




(9-32) 


the constitutive equations (9.5) through (9-7) yield (without further 
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approximat ion) 


Te 0 * 2v) + 3N 


(9.33) 


In our calculations the stress at the onset of necking was always 
at a level of stretch at which (i) (9-32) was a very good approximation, 
and (ii) the influence of v on h\ j“ was negligible. 

The dimensionless stress (r^Ay*) arises naturally in the 
analysis of Hill and Hutchinson. When the material (9.5)-(9-7) is 
incompressible, considering (9-9), we find 


7T Nt 


33 \ N 


m 


- *N In (X/X ) 


(9-3M* 


where X* (aV2a 3 ) is the stubbiness-j and X = (aV2a 3 ). It is well 

o o o 

known that the maximum lead occurs when the tangent modulus (as defined 

by . 39 or 9- 31 ) falls to equal the stress ; (i .e. when (t 33 /^*) = 1 ) . ** 

No bifurcation can occur before this point [23]. 

The tangent modulus continues to decline after the maximum load, 

so 4y : *<T 33 , H « y; that is, the tangent modulus is much smaller 

rnd x I o 3 0 

than the shear modulus in the neighborhood of the bifurcations points. 
As such, the critical stress (T 33 My") may be found by the asymptotic 
formula*** 


it 

note the latter equation is independent of any elasticity. 

& J* 1 

this can be shown by computing P*0 in equation (A.3M. 
see Hill and Hutchinson [59], equation 6.8 therein. 
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(T 33 /4y*) * i + 


sin 2y 


(9.35) 



where y*nnrX, m an integer. A$ (2y*/y)->-0, this formula reduces to 
that of Cowper and Onat [61] for a rigid-plastic solid. We could use 
(9*3*0 to eliminate either (T 33 /4y*) or X from (9*35) to get eigenvalue 
equations for X or (t J 'V I ty'), respectively, but for clarity it is bet- 
ter to plot (9.3*0 and (9-35) independently In the X- (x 33 Ay*) plane. 

The critical configurations in plane extension are then identified as the 
points at which the curves intersect. This is the approach we take, 
marking the critical configurations found by application of the finite 
element algorithm on the same plot. 

For numerical study we take Young's modulus E ■ 6.895 x 10 MPa, 

v-(l/3), and Ty-3M.75 MPa. Six individual cases are considered, 

corresponding to values of the hardening exponent N»A and N*8, for 

initial slendernesses of 1/X «2, 3* and *». These same six cases were 

o 

studied by Burke and Nix [60]. The problem may be treated in two 
parts: (i) generation of the solution for homogeneous extension, and 
(ii) location of the critical configurations through which the specimen 
passes in the course of homogeneous extension. 

To generate the solution for homogeneous extension for all six 
cases it is only necessary to find the solutions for extensions of a 


174 


unit cube of the two materials involved. These solutions consist of 
a sequence of configurations through which the specimen passes in the 
course of plane extension. It is anticipated that the bifurcation 
analysis will be sensitive to small variations of the stress and stub- 
biness, so an accurate integration of the homogeneous extension is 
essential. We use one element. For the material whose hardening 
exponent N equals four, a single RK2 step brings the material from the 
stress free state to the yield surface. This is followed by 30 RK*4 
steps to bring the specimen out to nominal stretch 1 = 1.0*4. Subsequent 
steps are (all RK*4) stretch increments Al= 0.002 out to 1 = 1.10, fol- 
lowed by Al = 0.01 out to 1=1.145. For the material whose hardening 
exponent N is eight, *40 steps were taken from the yield surface to get 
out to 1 = 1.05, followed by stretch increments A 1 =0-01 out to 
1 = 1.25 (all RKk) . For both N = *4 and N=8 post maximum load calcula- 
tions were repeated using Al= 0.002. This gave not only an accuracy 
check, but also a refined sequence of configurations in the neighbor- 
hood of the bifurcations points. One result of this increase in 
accuracy (over the example among the homogeneous deformations) is that 
the stress, as a function of the nominal stretch, is found to be 
virtually indistinguishable from the stress in the incompressible 
material: onlv the thickness of the specimen varied appreciably, .and 
this difference is not distinguishable on the figures which follow. 

In the second part of the problem, location of the critical 
conf igurgt inn* , estimates for thos« locations were provided by the 
results of Burke and Nix [60]. The global stiffness matrix [K*] 
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was formed using a uniform mesh of eight-noded elements, three elements 
through the (half-) thickness a' and twelve elements along the length 
2a^. Then the first eigenvalue (smallest in absolute value) was 
found.* The first eigenvalue has the physical significance of force 
required to produce 'unit' necking. The search for the critical con- 
figuration was extended towards or away from the maximum load point as 
this eigenvalue was negative or positive. The critical configurations 
are given In Table I. Linear interpolation was used to get configura- 
tions intermediate to those found in part one of this problem. In 
Figures 27 and 28 the present results are compared to the closed form 
results of Hill and Hutchinson [59] and to the numerical results of 
Burke and Nix [60], for N»k and N*8, respectively.** 

In each of Figures 27 and 28 the three nearly vertical lines 
represent the closed form solutions (9-3k) fo r X q = (1/2), (1/3), and 
(1/k). The stubbiness X decreases and the stress (t^/ky*) rises as 
extension progresses. Along the loading paths bifurcation is first 
possible when the curve (9-35) (for m*l) is encountered.*** Special 


this. 


* 


IMSL (FORTRAN) Library subroutines were 


used to accomplish 


Burke and Nix [60 ] do not reportjthe critical configurations, 
so it was necessary to reconstruct from information they do give. The 
stubbiness X can be recovered from equations k6a, k6b, and Table 2 (in 
that paper). The stress (T33/kp*) is then assigned according to 
(9.3k) — in this paper. We remark that their results are considerably 
more accurate than entries in (their) Table 2 would have the reader 
believe; the error In that table results f rom mi s-appl i cat ion of the 
asymptotic formula (7-6) from the paper of Hill and Hutchinson [ 59 ] - 


for m 


to show the proximity of the next bifurcation point, (9-35) 
2 has been plotted also. 
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Figure 27. Results of Necking Ana1ysis--N = A 



STRES 



symbols have been placed on the plots to Indicate the critical 
configurations found by the present finite element method and by the 
velocity-based finite element method used by Burke and Nix [60], 

Perhaps the most striking feature of these figures is that In 
every case the present method indicates bifurcation at a stress below 
the critical stress for an Incompressible body, while the results of 
Burke and Nix indicate the opposite. While the velocity-based method 
is certainly providing an upper bound for the critical stress, it would 
be Imprudent to assume from this single result that the stress-based 
method is leading to a lower bound. In fact, counterexamples to such 
a supposition are abundant. As an experiment, the Poisson ratio was 
varied between v*0.1 and v =* 0 . ^9 . It was found that as v increased, 
bifurcation was delayed. This is certainly not proof, but it supports 
the idea that the bifurcation stress of the incompressible body is an 
upper bound for the bifurcation stress of compressible bodies. In 
that light, the numerical results found by the present method must be 
more accurate than Figures 27 and 28 indicate, and the question of 


whether they lie above or below the correct critical stress (for 
v «y) is still open. Moreover, the present result was obtained using 
only 36 elements (239 unconstrained velocity parameters), whereas 75 
elements (1*79 unconstrained velocity parameters) were required In 


the velocity-based analysis to obtain a result less accurate.* 


Since a precise account of the critical configurations is not 
given by Burke and Nix {60], It is impossible to say whether or not 
the difference In the results can be attributed to the accuracy with 
which the extension was Integrated. 
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One could ask whether or not the symmetry we Imposed on the 
velocity field (see Figure 26} resulted in the suppression of an anti- 
symmetric bifurcation mode, such as one that formation of a shear band 
might give rise to. This Is unlikely, since the bifurcations all occur 
at stress levels too low to support the formation of a shear band in 
the incompressible material. That stress, for incompressible materials, 
is (see [ 60 ] ) 

(t 33 Au*) > V^TTT (9.36) 

Finally we note that the present bifurcation problems are well 
conditioned in the senses that (?) no other bifurcation points are 
close by, and (ii) the loading path and locus of critical stresses 
cross at a large angle. The results of Hill and Hutchinson [59] pro- 
vide a basis for evaluation of the performance of a finite element 
algorithm under much more demanding circumstances. 

A Parameter Study . In this section we investigate the sensi- 
tivity of the results of the previous bifurcation analysis to 
variations of number, shape, and type of elements in the finite element 
mesh. The 'type 1 of element is determined by the number of boundary 
nodes, the number of spin parameters, and the number of stress rate 
parameters. It Is desired that a 'small' change in the character of 
the mesh result in a 'small' (or at least predictable) change in the 
approximate solution. Moreover, it Is hoped that the accuracy of that 
result will increase monotonical ly with the computational effort, 
measured roughly by the number of unconstrained velocity parameters. 
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Finally, it is very important in practice to know whether an 
approximate load is an upper or lower bound for the actual load. 

Four different types of elements are considered. The first two 
are eight and four-noded elements satisfying all of the rank conditions 
given in Chapter VII. Third, a four-noded element whose spin field is 
of polynomial degree one less than the stress-rate field is considered. 
Finally, a constant st ress--constant spin element is considered. In 
all of the cases that follow, the configuration and stress of the body 
are fixed at values given in Table 2. The material is identical to 
that considered in the previous section with hardening parameter N=*8. 
Though the precise eigenvalue of the configuration (Table 2) is not 
known, it is believed to be small and positive. 

In Figure 29 the smallest eigenvalue of the global stiffness 
matrix [K' r ] is plotted as a function of the total number of uncon- 
strained velocity parameters. The finite element meshes were made up 
of uniform eight-noded quadrilaterals with 21 stress parameters and 
6 spin parameters. Complete data is given in Table 3- As can be seen 
in the Figure, the eigenvalue is quite insensitive to the particular 
arrangement of elements, depending almost exclusively upon the total 
number of degrees of freedom In the mesh. Every estimate for the 
eigenvalue was positive. 

In Figure 30 the smallest eigenvalue of the global stiffness 
matrix [K*] is plotted as a function of the total number of uncon- 
strained velocity parameters. The finite element meshes were made up 
of uniform four-noded quadrilaterals with 13 stress parameters and 3 
spin parameters. Complete data is given In Table In contrast to 
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Figure 30. Necking Eigenvalue for Various Finite Element Meshes--Four Node Elements (NT=13, NW=3) 


the eight-noded elements, the four-noded elements show a marked 
sensitivity to mesh arrangement. They are also substantially stiffer. 
In spite of the existence of preferred meshes, the eigenvalue estimate 
appears to improve as the degrees of freedom in the mesh increases.* 
Every estimate for the eigenvalue was positive. 

In Figure 31 the smallest eigenvalue of the global stiffness 
matrix [K' r ] is plotted as a function of the total number of uncon- 
strained velocity parameters. The finite element meshes were made up 
of uniform four-noded quadrilaterals with 13 stress parameters and 1 
spin parameter. Complete data Is given in Table 5. The stress-rate 
field on this element contained linear terms, while the spin field was 
a constant. Thus, angular momentum balance is generally satisifed 
only in the mean by this element. The reader should note the dramatic 
increase in stiffness of this element, as well as the relative insensi- 
tivity to mesh arrangement. Every estimate for the eigenvalue was 
pos i ti ve. 

From these three examples it is readily seen that the compliance 
increases with the number of kinematic degrees of freedom (velocity 
and spin parameters), just as it would in a velocity-based finite 
element algorithm. However the type of element makes a bigger dif- 
ference than does the number of elements in the mesh. This is 
important from the point of view of efficiency, since it means that 


the reader should note that the angularity of the lines con 
necting meshes in a sequence would be reduced if a refined sequence 
were used in the plot. 
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Figure 31. Necking Eigenvalue for Various Finite Element Meshes--Four Node Elements (NT* 13. NW = 1 ) 





it is better to go to a 'higher order 1 element than to refine the 
finite element mesh to achieve greater accuracy. It is also clear from 
these examples that the finite element algorithm is not providing a 
lower bound for the eigenvalue.* 

We would expect that decreasing the number of stress parameters 
would have just the opposite effect that decreasing the number of 
kinematic parameters had--decrease the stiffnesses (see Plan [ 3 ]). 
However, we found no difference in the necking eigenvalue for four- 
noded elements when 13 or 21 stress parameters were used, for either 
1, 3, **» or 6 spin parameters. This result is summarized in Table 6. 

As a final example we consider meshes of four-noded elements, 
each with 5 stress-rate parameters and 1 spin parameter. Tach 
element has two kinematic modes, but when the global stiffness matrix 
[K*] is assembled and the kinematic boundary conditions enforced, 
these modes disappear. The element is interesting because the 
equations of compatibility and angular momentum balance are satis! fed 
precisely on the interior of each element. As such, a velocity field 
of the form 

y * (a' x + b'\)e. (9-37) 

may be integrated on the interior of each element. However such a 
velocity field is incapable of matching the boundary velocities of a 


as might have been inferred from the results of the 
bifurcation analysis. 


G-3 
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four-noded quadrilateral except in the mean sense (since they contain 
an 'xy' term). Therefore, across element boundaries this velocity 
field must generally be discontinuous. 

According to Washizu,* one of Rayleigh’s principles states that 
"if the prescribed boundary conditions are partly relaxed, all the 
eigenvalues decrease." In using the elements above, the actual problem 
has been replaced by a problem in which interelement velocity con- 
tinuity has been partly relaxed. Thus, it is not surprising that for 
some meshes the finite element algorithm overestimates the compliance 
of the body. However, as the element mesh is refined, it can be 
imagined that the disparity between the boundary velocity and interior 
velocity diminishes (how rapidly would depend upon the particular 
sequence of meshes). In Table 7 two sequences of meshes are detailed. 
Any smooth velocity field could be approximated to any degree of 
accuracy by continuation of either of these sequences. In Figure 32 
the sequences of eigenvalues corresponding to these two sequences are 
plotted. These simple elements apparently are converging to the same 
value as all the other elements at a rate matched only by the eight- 
noded 'high-order' elements. But the most striking feature is that 
one of the sequences of approximate eigenvalues is converging from 
above while the ©♦'her from below. 

The natural tendency would be to attribute this behavior to 
the presence of kinematic modes on the element level. However, for 
sufficiently distorted meshes of other-wise well-behaved elements, 

Washizu [36], p- J*8. 
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This example supports the idea 


similar behavior can be observed.* 
that the present method does not necessarily lead to an upper or lower 
bound. 

Growth of a Void in a Viscoplastic Medium 

In this final example we examine the growth of a void in a 
hypoe last i c/vi scopl ast i c medium. This problem has been studied 
(numerically) by Burke and Nix [62], who treated the material as 
rigid/viscoplastic. We present the problem as a further demonstration 
of the performance of the finite element algorithm. The material 
exhibits stress relaxation, so the forward gradient scheme must be 
used to stabilize the time integration. Thus, only a qualitative 
picture of the stress and deformation can be expected of our analysis. 
Nevertheless, the present results agree quite closely with those of 
Burke and Nix [62). 

The motion is assumed to be plane strain, and throughout the 
body is a doubly periodic array of cylinderical voids. Due to the 
symmetry we need analyze only one quadrant of one rectangular cell of 
the body. The finite element mesh and boundary conditions are 
described in Figure 33. 

Burke and Nix motivate their study by explaining that certain 
theories for the initiation of creep fracture suppose that the growth 


For a (1 x 10) mesh of eight-noded elements an eigenvalue of 
-30.3368 was found; it was not determined whether this value was 
simply erratic (due to the irregular element shape), or whether 
similar negative values could be found for 'nearby' meshes (such as 
1x9 or lxll). 
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Figure 33. Finite Element Mesh and Boundary Conditions for Void 
Growth Problem 
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of voids can be "attributed to the inhomogeneous plastic deformation 
of the surrounding grains." Furthermore, "finite fracture strains can 
be predicted only when a void lies in the neighborhood of another 
void." Such a study necessarily involves a number of special cases. 
For our purposes, that of demonstration, only one case is taken. 

The problem has been analyzed in three parts, much as the pipe- 
creep problem was. in the first part the cell is brought rapidly from 
the virgin state (s tress- free) to a state of purely elastic strain. 
This is accomplished by a single RK2 step. In the second part, 
relatively small time steps are taken while the stress relaxes from 
the elastic distribution to a nearly steady creep distribution. In 
the third part, time steps are taken which produce 1% nominal 
elongation of the cell in each step. To stabilize time integration 
in the second and third parts the forward gradient scheme is used, 
the stability parameter 0 set as 0=1/2 and 3A> respectively. Con- 
sistent with our earlier discussions regarding use of the forward 
gradient scheme, only the Euler time stepping scheme has been used in 
the second and third parts of the problem. 

The material model is identical to (9-11): 
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This model corresponds to that of Burke and Nix [62] with (their) creep 
exponent n* 1. The fluidity y is set as y« ! x 10 (psi-sec) . The 

velocity at the top of the cell (see Figure 33) was adjusted so that a 
specimen with no void would experience a homogeneous constant stretching 
e 1 ^ of e 1 1 * 0.25 x 10 ^ sec ^ . Since the material was treated as rigid/ 
viscoplastic in [62], our choice of elastic constants is somewhat 
arbitrary. We have taken Young's modulus E»3xl0^ psi and Poisson 
ratio * 0 . 4, so the material is like mild steel in its elastic response . 

In Figures 34, 35, and 36 the contours of stress T 1 \ mean stress, 
and stress have been plotted for L (the elongation of the cell) 

L 1.01. The stress concentration where the hole edge crosses the x^ 
axis is approximately 2.7.* This is quite reasonable since the 
theoretical value for an isolated void in a purely elastic medium is 
3.0 [63]. in Reference 62 an approximate value of 2.66 was found for 
the rigid plastic material. in Figure 37 the contours of effective 
strain rate 
compares very well to Figure 7 in [62]. 

in Figure 38 the deformation is traced from 1*1.0 to L=1.5. 

These deformations are physically tenable. We remark that no Indication 
of any numerical instability was observed in the course of integrating 
this deformation. 

In Figures 39, 40, and 41 the contours of stress T 11 , mean 
stress, and have been plotted for L*1.50. They compare very well 


A stress concentration of approximately 2.59 was observed for 
the elastically stressed medium. 


V 


2 P 

3 ~ 


E are plotted for L ** 1.01. Qualitatively this 
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to the stresses found In [62] (see Figure 8 there). We note that the 
stress concentration has dropped to 1.71. The stress concentration 
depends strongly on the geometry of the specimen; as such. It was 
observed to decline steadily throughout the deformation. In Figure ^2 
the contours of effective strain rate are plotted for L * 1.5. Again, 
the qualitative agreement with the results of Burke and Nix [62] is 
noted (see Figure 9 there). 

The present calculation was terminated at L= 1.5 because of the 
unstable traction boundary condition* at x^*0.0 and the edge of the 
hole, and the general breakdown of (total) traction reciprocity con- 
ditions on the interior of the cell. This problem is easily avoided 
by incorporation of traction residuals. 

We conclude by noting that in the present analysis only 56 four 
noded elements were used, as compared to 56 eight noded elements used 
in the analysis of Burke and Nix. Considering the agreement between 
their results and our own, the present method appears to have performed 
very well, in spite of the large disparity in the degrees of freedom 
of the finite element mesh. 


see discussion and footnote accompanying equation (3.37). 
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•CHAPTER X 


CONCLUSIONS AND RECOMMENDATIONS 
Concl us i ons 

In this work a new hybrid stress finite element algorithm, 
suitable for analyses of large quasistatic deformations of inelastic 
solids, has been presented. The feasibility and performance of the 
algorithm has been demonstrated in a number of example problems. 

Principal variables in the formulation are the nominal stress 
rate and spin. As such, consistent reformulation of the constitutive 
equation is necessary. This is discussed at length, as are alterna- 
tives to direct numerical inversion of constitutive matrices involved 
in that reformulation. 

The principal variables in most finite element algorithms for 
solids are either displacement increments or stress increments (as 
opposed to actual time derivatives). In problems involving elastic 
bodies the accumulated error of such an algorithm may be kept small 
by 'residual load' iterations; however, in problems involving inelastic 
bodies, the accumulated error of the incremental approach is like that 
of the Euler scheme for integration of ordinary differential equations. 
In the present work the notion of 'increments' has been discarded 
entirely. As a consequence, the finite element equations give rise 
to an initial value problem (which may be treated independently). 
Integration has been accomplished by Euler and Runge-Kutta schemes, 
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and the superior accuracy of the higher order schemes Is noted. 

It has been shown that there is an ambiguity Inherent in finite 
element methods based on complementary work and energy principles sur- 
rounding the appropriate definition of the velocity (or displacement) 
on the interior of the element. Heuristic arguments have been given as 
justification for the method by which those velocities were found in 
the present work. Those arguments indicate that mathematical con- 
sistency requires that special techniques be used to find the velocity 
on the interior of 'high order' elements, but that the velocity on the 
interior of a 'low order' element may be found by interpolation of the 
boundary velocities. 

In the course of integration of the stress (in time) it has been 
demonstrated that classical schemes such as Euler's and Runge-Kutta may 
lead to strong frame dependence. The problem can be traced to the 
integration schemes themselves. As a remedy, modified Integration 
schemes have been proposed. The potential of the new schemes for 
suppressing frame dependence of numerically integrated stress Is 
demonstrated by an example. 

Time integration of the stress in materials which exhibit stress 
relaxation is complicated by the necessity that one take very small 
time steps In order to avoid numerical instability. The applicabi lity 
of explicit and implicit forward gradient schemes to Improve stability 
of Integration in large deformation problems has been Investigated. 

These schemes are known to be both stable and accurate in problems 
involving small deformations. It has been found that in large deforma- 
tion problems the schemes are indeed stable, but potentially inaccurate. 
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The capability of the stress-based finite element algorithm for 
extremely accurate bifurcation analysis was demonstrated. Moreover, it 
was shown that one could expect the result of such an analysis to be 
stable with respect to variations of the finite element mesh, so long 
as the same type of element was used in every mesh. If the type of 
element was varied, the result changed in a (qualitatively) predictable 
manner. It was demonstrated that the method did not necessarily lead 
to an upper or lower bound for the critical load. 

Finally, it was made evident through examples that stresses 
obtained by the present method were of exceptional accuracy; much more 
than could be expected of a velocity-based algorithm. Traction 
boundary conditions and the traction reciprocity conditions were met 
with a correspondingly high accuracy, though their accuracy could have 
been improved by incorporation of residuals (to keep the accumulated 
error smal 1 ) . 

Recommendat ions 

The principal defect of the algorithm presented in this work is 
its inability to accurately integrate the stress in bodies which exhibit 
stress relaxation, unless of course, the time steps are kept inordi- 
nately small. Evidently, this defect is present in the algorithm of 
Kanchi et al [^9], though no mention is made of it. This appears to be 
the only other application of a 'forward gradient' technique to finite 
deformation problems. In the context of small deformation problems, 
the error of (generalized) gradient techniques was studied numerically 
by Argyris et al . [19]. On that basis they concluded that 'midstep 
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weighting techniques' (that is, 0*i) were accurate. However, the 
present results Indicate that the accuracy of generalized gradient 
techniques is highly problem dependent, and that no one choice of 9 
assures optimal accuracy. A minimal requirement to be made of any scheme 
for stabilization of numerical integration is that it give correct 
results under steady, or nearly steady conditions; this requirement is 
not met by the generalized gradient techniques currently available. It 
is recommended that an effort be made to develop schemes for stabiliza- 
tion of numerical integration of stress for finite e lement- i n i t ia 1 
value problems whose accuracy can be proved. 

Secondly, it was demonstrated that the present method could not 
be relied upon to give either an upper or lower bound for the critical 
load in a bifurcation analysis. It is suspected that the character of 
the approximate load obtained by the present method may be linked to 
the rank conditions (7.69) and (7-73). Further research, both from the 
mathematical point of view and numerical point of view is needed before 
a practically applicable criteria for critical load characterization 
can be given. A mathematically accurate discussion of the problem of 
assigning the velocity on the interior of an element (for hybrid stress 
methods) might aid in this characterization. 

Finally, the materials considered in this work all have the com- 
mon property that their constitutive equations are isotropic functions. 
Though the generalization of the present method to anisotropic material 
behavior is conceptually straightforward, there may be special problems 
in an implementation. The performance and special problems of the pres- 
ent method, when applied to anisotropic materials, should be investigated. 
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APPENDIX A 


DIRECT, DYADIC, AND INDEX NOTATIONS FOR TENSORS 

Let r denote the position vector in We write X* for an 
arbitrary (smooth) system of coordinates in E. Then a triad of vec- 
tors e | , called the 'natural base vectors' of the coordinate system, 
are defined by the equation 


S| ■ — 

ax' 


(A.l) 


We assume that these vectors are linearly independent; that is, 


(£j x e^) • e. / 0 


(A. 2) 


Then a conjugate triad of base vectors Is defined by the equation 


s' * «J 


1 if I = J 


0 otherwise 


(A. 3) 


Any vector in £ may be represented as a linear combination of 
the base vecturs e ^ or e* : 


I I 

v = v e . * v ,e 


(A.M 


The cor-iyonents v* ar c called ' cont ravari ant 1 and the components v ( are 
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called ‘covariant. 1 They are defined by 


I 

v 



S| 


(A. 5) 


1 2 3 

To every triplet <v ,v ,v > there corresponds a unique vector in E, 

and vice versa, for any linearly independent set of base vectors ej. 

This important fact is the basis of the index notation for vectors 

and tensors. If a linearly independent set of base vectors is given 

(for example by specifying a coordinate system), then no ambiguity 

12 3 I 

arises in writing <v ,v ,v >, or simply v , for v. In continuum 
mechanics one frequently is forced to work with more than one 
coordinate system at a time; then the dyadic notation, indicated by 
(A. A), is more convenient. 

The relation between the contravariant and covariant components 
is found by 'dotting 1 the representation (A. A): 


v 1 » (e 1 , e J )v J ; » (e ( • Sj)v J . (A. 6) 

Just as the triad «t| constitutes a vector basis in E, the dyads 
are a basis for second order tensors in E. Any second order 
tensor may be represented as a linear combination of these dyads, for 
examp 1 e 


_ . _| J _ I J / I J\ T J , I . 

T ■ T e j e j ; T ■ (e e ) ; T « e • (e • T) 


(A. 7) 


Other components for T may be defined for the dyads 
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and (e ( e^), and they are called 'covariant,' 'mixed,' and 'mixed,' 
respectively. To every 3x3 matrix there corresponds a unique second 
order tensor In E, and vice versa, for any linearly independent set of 
dyads. Therefore the index notation may be used for second order 
tensors also. . r 

A fundamental tensor In E is the Identity tensor, defined as 
that tensor with the property - 


v * I • v 


(A. 8) 


for every vector in E. We may represent I as 


-I J 
I » 5 j£|S , 


where 5 ^ is the matrix of components of I for the dyad (eje^). Accord- 
ing to (A. 5) we may represent v by v*e|, so (A. 8) may be set in the 
form 


(v 1 - 6jV J )e , - 0 . 


(A. 9) 


For (A. 9) to be satisfied for arbitrary v 1 it is necessary that 6 1 j be 
defined as 


6 V 


1 if I - J 


0 otherwise 


• (e • e j ) 


(A. 10) 
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Components of I may be obtained for the other dyads in a similar manner 
as 


6 IJ « (e 1 • e J ) ; - (e t • e^) ; 5j J - (e ( • e J ) . 

(A. 11) 

In some applications I Is called the ’metric tensor. 1 

The relation (A. 6) between cont ravar i ant and covariant vector 
components may be written as 


I 

v 




(A. 12) 


Similarly, the relations between the various components of a second 
order tensor may be written as 


T 1 J - S KJ T ! k - $ K 't k J - 5 KI 5 U T KL • 


(A. 13) 


The apparent rule for raising and lowering indices may be shown to be 
valid for tensors of all orders. 

When differentiating vectors along a coordinate line, one must 
take into account not only the rate of change of the components, but 
also the rate of change of the base vectors themselves: 


3 / \ 3 i K 3 / \ 

— r (y) “ — r ' v > e u * v — r • 

3X 3X 3X 1 K 


(A.UO 


The derivative of a base vector Is a vector itself, and may be 
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represented as a linear combination of the base vectors: 



y ik?j 


(A. 15) 


The components In the representation (A. 15) are called Christoffel 
symbols (of the second kind); they are defined In the same manner as 
components of any other vector, by use of (A. 5). Using the representa- 
tion (A. 15), the derivative of the vector in (A.lA) may be written 




+ v 


K J \ 

y ki ; 


(A. 16) 


The coefficient in (A. 16) Is called the 'covariant derivative' of the 
component v J . Covariant derivatives of the components of second order 
tensors may be defined in' a similar manner. It should be noted that 
the Christoffel symbols are not components of third order tensors. 

The operators GRAD, DIV, and CURL may be represented as the 
vector operators 

GRAD v - Vv 


D l V v = v • v 


(A. 17) 


CURL v - V x v 


where 7 is the symbolic gradient operator: 
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(A. 18) 


V 



In dyad notation 


3 



(A. 17) may be written out as 


GRAD y - Vy - (v J ), , 5 '^ 

DIV y « V • y - (v J ) , (A. 19) 

CURL y * V x y » Cv l ), J e J |(< e K 

where ( ),| denotes the covariant derivative with respect to X*, and 
e'j^is the alternating tensor, defined by 

e * ■ (e * x e ) * e *5^(e x e ) • e = 6^e (a 20 ) 

JK -y -K l -N -K 0 NJK * 


A convenient summary of formulas of vector analysis is given by 
Spiegel [50]; an extensive treatment of the subject is given by 
Phillips [51]. The remainder of this appendix is devoted to special 
notations used in this work. 

The special notations used in this work are summarized in the 
formulas below. In accordance with (A. 7) we write a second order 
tensor T as 


T t 1 *1 

I “ T £|£j 


T 'j£|S J 


T J I 

T l S Sj 


T 1 J 

T l J- S 


The transpose Is given by 
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(A. 21) 


T T T U yl J - J I x J 

I " T SjS| - T J* S| - T , Sj£ ■ T ,jf S| 


A fourth order tensor E may be written out as 


r m C IJKL 

E E • 


The scalar product of two second order tensors S and T is 


S : T * S ,J T 


IJ 


The product of a fourth order tensor E and a second order tensor 


J'l- eIJKIt kl*i*j ; 


T , T .KLIJ 

I 1 J ■ T KL E S,£j - 


The product of two fourth order tensors D and E is 


_ _ -KLMN I J 

5 ;E -°IJKL E S S Sh'n • 


Finally, differentiation by a tensor is defined 


3R 3R 

3T “ 3T (J -l s J ; 


3S 

W ■ 3T kl - I - J- K-L 


(A. 22) 

(A. 23) 
T is 

«S» 

(A. 24) 

(A. 25) 
(A. 26) 
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ALTERNATIVES TO DIRECT NUMERICAL INVERSION 
OF THE CONSTITUTIVE EQUATION 

Formulation of the stiffness matrix involves inversion of the 
(9x9) constitutive matrix W at each quadrature point in the body, each 
time the stiffness matrix is evaluated. In practice it is found that 
these inversions figure significantly in the total computational 
effort. In this appendix we investigate the possibility of (1) analytic 
inversion of matrices of the form (4.21), and (2) approximation of the 
inverse of a matrix W = V-T when V ' and T are known. The reader is 

** as at M as 

referred to the articles of Rivlin and Ericksen [25], and Rivlin [52] 
for discussions of representations of symmetric isotropic matrix 
functions . 


Analytic Inversion of the Constitutive Equation 
We begin by inverting the counterpart of (4.21) which arises in 
plane stress and plane strain problems. We consider a symmetric 
isotropic matrix function of the form 


a* - V : e + I 

where a*, £ , and E are (2x2) matrices, and 
2 

X + 2l, i 

I , J*1 


(B.l) 


(B.2) 
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t 


£ - > n'z, 

i - 1 


(B.3) 


5i " 1 1 


52 " I- 


(B.M 


If V is invertible, then we may write (B.I) as 


-1 •* 1 

e * V : a* + Z 


(B.5) 


where 


t 


I, J-l 


A . z . + 2M I 

— | — J 


(B.6) 


and 


Z 1 * -v " 1 : Z . 


(B . 7) 


Of course V and V ^ satisfy 

* 


V : V _1 - V _1 : V - ! 


(B.8) 


so multiplication of (B.2) and (B.6) must yield 


( , IJ.KL, 

x A ( !j : !k 


L J? 1?L 


) + 2MX IL + 2yA IL |z,z, + Hum] I . 1 . (B.9) 

5* s 


(in which the summation convention is used). From the linear 
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Independence of the basis z.z., I we conclude that 

~ I ~L sj 


= 1 ; (B. 10) 

X U A KL (z j : z K ) + 2MX ,L + 2pA ,L * 0 . (B.ll) 


From (B.10) we get 


2M = 1 / (2y ) . 


(B . 12) 


Elimination of 2M for. (B.ll) leads to a (2x2) matrix equation for A IJ : 
[[X][Z] + 2y[ I ]] [AJ - - ± [X) (B. 13) 


where 



> 

x ,2 ‘ 


- 

— 1 

> — 

> 

to 

M - 

■ X 2 ' 

X 22 . 

f 

IA] - 

.A 21 A 22 J 


and 


I J 


(z 


Zl 




The solution of (5.13) is easily found by Kramer's rule: 
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a2 ' ” 2^A [^"X 22 - A ,2 X 2, )(trT)I 

^.^[ u "^.x'v T )( trI ,]-4i 

A - (A n X 22 - A 12 A 21 )(2(trj 2 ) - (trr) 2 ) 

+ 2y(2A U + (A 12 + X 21 )(trr) + A 22 (trr 2 )) + Ay 2 . 

-x ««# 


Thus 2M and A*^ in the representation of V ' (B.6) are all explicitly 
determined by y, A IJ , and the stress. 

The problem of inverting V (A. 21) for general problems may be 


attacked in precisely the same way as for two dimensional problems. 

- 1 ” 

For V and V we write 


V 

* 


x hh 


2y $J 



a'^z.z . + 2M 1 4> 

— 1 — J as 


(B. 15) 


where 


5j 


52 


53 
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and 


[0. ] . .. . ■ 6. , 5 . 

1 1 J i jk 1 i k 1 j 


'Vljkl ' i(T i'k 5 lJ + S lk T ij> 


'♦jlljkl ' J(s ik s ij +s ik s lj> • 


The equation V ^ : V «* I gives the following relation among the X 1 ^, 

I A I J i J 
y , A , and M : 


[x IJ A KL (z l : z,) 5k 2 , f 2y M A KL z K (z L :$„) 


(B.16) 




It is tedious but straightforward to resolve the expressions (z. : 0 .) , 

(0, : z,), and (0. : 0.) in the basis z.z., 0,. The formulas given by 
« I % J -j; I 

Rivlin [52] (generalizations of the Hami 1 ton-Cay ley theorem) are 
particularly useful. 

Resol utions : ( B . 1 7 ) 


*i ! £i *£rji *£i ' " ,l2 ’ 3 • 

*2 : !i ’ £i : *2 ’ £2 



£2 



£3 
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$2 : £3 " *3 : I2 * ^ det I ')*J + *( tr Vh 

ti : ii ■ Jl : b * 53 

*3 : i2 ’ ?2 : *3 ' *2 : 53 

$3 : Z3 ■ - (det x')z 2 + i(tr s)z ^ 

!l : " $| : $i " $| I - 1 ,2,3 

*2 : *2 ’ - J +3 ' i(tr 5*5i2i + *h52 

+ 1(ZjZj + Z^Zj) + 4(tr s)^ 

<j>, : - <j>, : 4)- - i(tr s)<j>, + 4(det T* ) I 

+ 4(z^z 2 + z 2 z 3 ) + £(det T 1 )zjZj 
4> 7 : " r ( tr sH* “ i(det t') 4>_ - i(t r s) 2 i 

.-3 '~5 *♦ -» mi -» »Z ~ « 

* ^353 ' i(tr ! )( 535 i + 5i53> 

+ i(det + ZjJ,) + j(tr s) 2 ^,, , 

After resolution of the terms in (B . 16) into the basis of the 
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representation of V (i.e. ZjZj and 4> ^ , the coefficient of I Is set to 

one and the coefficients of the other terms are set to zero. The 

.1 ..I 


coefficients of the 4>. 

as * 

i nvol ve 

the generalized shear 

modul 

and 

.tress 

only: 




1 : 

*u' 

y 2 (trs ) + y 

^ (det t‘ 

) -2y 3 (trs) 2 +y 2 (det 

t ) 

^2 : 


ijy 1 + y^(trs) 

y 2 (trs) - 2y^(det x 

') 

*3 : 


-2y 2 


.^y 1 + 3y 3 (trs) 




In general (B . 18) could just as well be solved numerically as analyti- 
cally. In a special case of great practical importance though, when 

2 7 11 

y and both vanish, it follows immediately that 2M «*(l/2y ), 

2M =2M J = 0. In any case, the remaining equations form a (3x3) matrix 

equation for the A*"*. When 2M^ * (l/2y^) and the other M* vanish, the 

equation for the is of the same form as (B.13): 



(B. 19) 


where the matrices [X], [A], and [Z] are the identical counterparts of 
those in (B.13). 

Though an analytic expression for V ' was not found in the 

31 

general case, the numerical problem of inverting a (9x9) matrix was 
replaced by the problem of solving (at most) two (3x3) matrix equa- 
tions. Finally we remark that no assumption as to the symmetry of the 
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[X] matrix was made 


Approximation of W 

Suppose that we are given the constitutive equation for a body 
In the form 


e » V* 1 : o* + e P 


(B. 20) 


and we wish to obtain the form 


-1 • VP 

e - w : r + e . 


(B.21) 


Since a* «■ r + T : £ , we can get the impl i c i t equation 


« 1 . • . P 

e»V : (r + T : e) + e 


(B.22) 


directly. If we try to solve (B.22) for e (when r is assigned) by 
i te rat ion , 


e N+ ' ■ v"' : (r + T : E N ) * d 


(B.23) 


then we are led to define W.J as 


w" 1 - V H + v" 1 : T : v' 1 + . . . + (V’ 1 : T) H ; v’ 1 

>»N a 2 as as ~ 


(8.24) 


25 W 


• 1 N+l 

and the first neglected term is (V : T) : e, where (N + l) is the 

H Z " ** ' 


222 


I HIM II II 



number of terms in the series (B.2*0. If the eigenvalue of (V ' : T) 

a ss 

whose absolute value is greatest is of absolute value less than one, 
then the remainder vanishes as N + », so the series converges; that is 
W '-*W The eigenvalues of T coincide with those of the true stress, 
and the eigenvalue of V 1 whose absolute value Is greatest is a shear 
compliance. Thus, for metals in the elastic range W. can be expected 

as I 

to be. in error by less than 0.01%. We find that in practice if more 
than two terms in (B.2*0 are needed, it is more efficient to compute 
W ^ by some other means. The main appeal of (B.2*0 with two terms 

as 

taken is in large deformation — small strain analyses, such as in 
structures . 

Construction and Inversion of Constitutive 
Equations for Plane Problems 

We first indicate the class of problems which may be considered 
'planar.' The class consists of those problems in which the true 
stress t, the (general) stress rate s, and the stretching e are of the 
fo rms 


T a6 e e + t 33 e e 
T -a-B -3*3 


•aB_ _ *33 


i -a-B " 4 -3-3 


a3 33 

e * £ e e 0 + e^ J e,e, 
- ~ i 


(B.25) 


where the Greek indices range from 1 to 2. Substitution of (B.25) 
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Into the constitutive equation 


s ■ V : z + Z 


(B.26) 


yields the component equations 



r V aS y5 V aS 33 
Lv 33 y6 V 3 3 3 3 



(B. 27) 
(B.28) 


The tensors V" and I", defined as 


V" * V a6 y5e e e Y e 5 

a- “CTP* ” 


i" - z oS s a s s 


(B.29) 


are 


necessarily of the forms (B.2) and (B . 3 ) when V and I are of 
the forms (4.21) and (A. 22), respectively. We define s", z , and T as 


i" - s“ 6 s a s 8 ; 5 " * £<lf W- l" * : <S“S B » ( SaS 6 > • 


For plane ~t r 3in (B.27) may be written 


•H .,11 " , r »l 

s * V : z + I 


; 33 - V 33 : e" + I 33 


(B.30) 
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For plane stress (B.27) becomes 


• n 

s 




+ Z 


,33 £ 33 . v 33 
33 ~ 


+ I 


33 


(B . 3 1 ) 


It is apparent that if we can write V ' in the form (B.2) then 
analytic inversion is possible. To find the necessary coefficients we 
set the components so that 


(-a-g) ; ~ : ^-y-6^ “ ^-a-6^ ; : ^-y-6^ = 0 (B . 32) 

where V is written for (A. 21) and v" for (B.2). We let and y* be 

% m 

the coefficients in (I*. 21) and 1 3 and m be the coefficients in (B.2). 
Then (B.32) leads directly to 

i" -x" - q <x ,3 + x 3 ’) *, 2 x 33 - (pu 2 *pV) 

l' 2 -X ,2 + p(X ,3 + p 3 ) - ,X 32 - qpX 33 ♦ p 2 
I 2 ’ - X 2 ' ♦ P (P 3 ' ♦ p 3 > - qX 23 - qpX 33 ♦ p 2 
l 22 - X 22 * p(X 23 * X 32 ) * p 2 X 33 
2m ■ 2p' + pp 2 + (p 2 - 2q)p 3 
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i 2 

where p-T^ and q ■ £(p - T ^g T ^g) • The planar Inverse Is found by 
putting 1 for X , m for y, p for trT , and (p - 2q) for tr(t ) In 
(B.14). Though (B. 33) and (B . 14) are algebraically complicated, their 
effect Is to reduce the problem of construction and Inversion of 
for planar problems to about ten lines of ordinary FORTRAN. 
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APPENDIX C 


SHAPE FUNCTIONS FOR VELOCITY, STRESS 
RATE, AND SPIN 

Shape Functions for Plane Strain 



Velocity Shape Functions 


N. » N, .e. + N. .e. 

-i l,i-l 3,i-3 


Four Noded Element: 


’ i(1 +K 1 )0 +rm.) 1-1,2, 3,^. 

,0 l * 5, 6, 7, 8 



0 1 - 1 , 2 , 3 ,*» 

I - 5 , 6 , 7, 8 


1CI < 1 , hi <1, 

c 2 ■ 1. C3 * 1. Ci, - “1 

n ] - -l, n 2 - - i , n 3 - 1 , - 1 
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Eight Noded Element 


( H, 


N, . 


! - 1 , 2 ,., . ,8 

i -9,10 16 


N, . 

3,i 


/ 0 


l H i-S 


1 - 1 , 2, ... ,8 

I -9, ?0 16 


1 i(i - C 2 )(l +nn.) 

id +«,)(! V) 


i « 2,6 
i - 4,8 


4(1 + £■£.) (1 +nn- ) (CC. + rin. - D f - l ,3,5,7 


h * h 


-1; C 2 -C 6 - 0; E, 


3 ^5 


h 2 - n 3 


-i; \ m n 8 » 0 ; n e - ^ - n 7 - +i 


5 '6 7 


Shape Functions for Spin 


^i * QW 13,i -1*3 + ^31,1 -3-1 


where 

9U 13.I ’ C 1 

« U 3.,1 *-'l 

If NW-3 add the following shape functions 
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3U I3,2 ' xc 2 
5W !3.3 ‘ zc 3 
® W 3 1 . 2 ■ " xc 2 
^ W 31 ,3 " ~ zc 3 

If NW= 4 add the following shape functions 


QW 


13,4 * xyc 4 


Q.W 


31,4 " ” xyc 4 


If NW = 6 add the following shape functions 


qw I 3,5 ’ xxc 5 
5W 13,6 ' zzc 6 
qu 31,5 * " xxc 5 
5U 3I,6 ' " zzc 6 

The constants were used to improve the condition of [H], 


Stress Shape Functions 


QJl " QT 1 1 ,1 SlS] + 0 + ^ T 13 ,T S 1 S 3 


: * -22,1 -2-2 + ^ 31,1 *3 S 1 + 0 + I US 


33,1 -3-3 


- 1 

• t * 

5T 31,2 - " 


QT 22.3 ' ' 
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- -1 


QT 


13,4 


"W 1 

For NT" 13 add the following stress shape functions 


QT 


11,6 


QT 


31.6 


-z 


ttT 31,7 " 


QT, 


22,8 


QT 13,3 * •* 
aT 33.9 ‘ 2 
aT 33.10 ’ * 


aT n,n 

QT 13.12 


m z 


• - Z 


QT 22,13 ’ 2 

For NT- 21 add the following stress shape functions 

* T ll,l* * * 5xx 


QT 


31,14 


■xz 


QT 


31,15 


« -XX 


QT 


13,16 


" - . 5*x 


QT 


33,16 


■ xz 


QT 


33,17 


xx 
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xz 


QT 

QT 

Q.T 

QT 

QT 

QJ 


11,18 

3M8 

11.19 

13.20 

33.20 

13.21 


*.5zz 


zz 

-xz 

,5zz 

-zz 


1x1 Gauss quadrature for 'constant' shapes; 
3x3 Gauss quadrature for 'linear' shapes; 

A x 4 Gauss quadrature for 'quadratic' shapes. 


Shape Functions for Ax i symmetric Deformation 


1 

r = x ; 



-1 " -r ; -2 " s 0 ; - 3 * -z ’ 

Components N. N, and N. . are identical to plane strain shape 

1,1 2,1 i , I 

functions. Spin functions are identical except r replaces x, e^ 
replaces e^. 

Stress Rate Shape Functions 



QT n , I + 0 + QT 13,i *1-3 
+ 0 + QT 22,i -2-2 + 0 + liT 31 , i *3*1 + 0 + ^ T 33.1 S 3 ? 3 
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where 


aT n,i 

c/r . 

* T 31.2 

c/r . 

QT 13.3 

c/r . 

QT 33,A 

c/r . 

* T H,5 

1 . 

QT 22,5 

1 . 

QT 31 ,6 

1 . 

QT 1 3,7 

c . 

QT 33,7 

-cz/r . 

^33,8 

cz/r . 

QT 11,10 

• z . 

QT 22,I0 

» z . 

QT 13.11 

■ cz/r . 

QT 1 1,12 

• (r-c)/c . 

QT 22,12 

- (2r-c)/c . 

QT 31,13 

- (r-c)/c . 

QT 13,1^ 

- (r-c)/c . 

QT 33, 1* 

■ -z(2r-c)/rc 

QT 33J5 

■ (r-c)/c 


ORIGINAL bm E f S 

0F POOR (QUALITY i 
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where c Is chosen to Improve the condition of matrix [H]. (3x3) 
Gaussian quadrature was used on this element. The constant c was 
assigned as the value of r at the center quadrature point on each 
element. 
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APPENDIX D 


TABLES 


Table 1. Critical Configurations in Plane Extension 


N« 4 

2aW 
o o 

a 1 

c 

2a c 

N 

(M O 
1- 

T 33 

C 

§ 

2 

3.6191 

13.93**2 

584.79 

1173-34 


3 

3.7800 

20.0053 

564.15 

1132.54 


4 

3.8388 

26.2625 

556.25 

1116. 87 

- 

N = 8 





- 

2 

A. 071** 

12.3385 

322.63 

646.19 

r 

3 

**.2852 

17.5819 

311.39 

623.93 

5 

4 

*♦.3533 

23.07**1 

307.26 

615.77 

- 


Table 2. Configuration for Parameter Study 


E 

2, 3 /a' 

1 

a 

2a 3 

T 22 

T 33 


0 o 






2 

**.0752 

12.3271 

322. M* 

645-83 

= 
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Table 3. Data for Figure 29 (8 Noded Element, NT ■ 21 , NW ■ 6) 


Mesh 

Symbol 

Degrees of 
Freedom 

Ei genval ue 

1 x 2 

□ 

15 

0.5341 

1 x 3 

□ 

23 

0.3221 

1 x 4 

□ 

31 

0.2416 

1 x 6 

□ 

*♦7 

0.1505 

2x4 

A 

55 

0.1850 

2x6 

A 

83 

0.1259 

2x8 

A 

111 

0.0957 

2 x 10 

A 

139 

0.0772 

3 x 6 

O 

119 

0.0925 

3x8 

O 

159 

0.0703 

3 x 10 

O 

199 

0.0568 

3 x 12 

O 

239 

0.0476 

A x 9 

0 

233 

0.0495 


Table 4. 

Data for Figure 30 

(4 Noded Element, NT - 

13, NW - 3) 

Mesh 

Symbol 

Degrees of 
F reedom 

Ei genval ue 

2x6 

D 

29 

0.3732 

2x8 

□ 

39 

11.2843 

2 x 10 

□ 

49 

22.6500 

3x6 

A 

41 

27.1690 

3 x 9 

A 

62 

0.1777 

3 x 12 

a 

83 

2.7302 

3 x 15 


104 

5.4067 

4x8 

o 

71 

9.5265 

4 x 12 

o 

107 

0.1059 

4 x 20 

o 

179 

1.8974 

5 x 10 

0 

109 

4.1738 

5 x 15 

0 

164 

0.0706 

5 x 20 

0 

219 

0.4415 
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Table 5. 

Data for Figure 31 

(4 Noded Element, 

NT - 13, NW - 1) 

Mesh 

Symbol 

Degrees of 
Freedom 

Ei genval ue 

3x6 

□ 

41 

633-077 

3x9 

□ 

62 

316.414 

3 x 12 

□ 

83 

228.714 

3 x 15 

O 

104 

197.089 

3 x 18 

□ 

125 

185.583 

1* x 8 

A 

71 

237.617 

x 12 

A 

107 

117.908 

4 x 16 

A 

1*3 

80.444 

4 x 20 

A 

179 

68.967 

5x10 

O 

109 

107.518 

5 x 15 

O 

164 

49.623 

5 x 20 

o 

219 

35.078 


Table 6. 

Necking Ei genval ue-- (2 x 6 Mesh, 

4 Noded Elements) 


NT - 13 

NT - 21 

NW - 1 

1253.19 

1253.19 

NW - 3 

0.3732 

0.3732 

NW « 4 

0.3732 

0.3732 

NW * 6 

0.3732 

0.3732 
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Table 7. Data for Figure 32 (4 Noded Element, NT » 5, NW ■ 1) 


Mesh 

Symbol 

Degrees of 
Freedom 

El genval ue 

2x4 

o 

19 

1.0948 

3 x 6 

□ 

41 

0.4043 

4x8 

□ 

71 

0.2074 

5 x 10 

□ 

109 

0.1269 

6 x 12 

O 

155 

0.0862 

2x8 

A 

39 

-1.2258 

3 x 12 

A 

83 

-0.2204 

4 x 16 

A 

143 

-0.0426 

5 x 20 

A 

219 

0.0003 
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Corrections to 


"Analysis of Large Quasistatic Deformations of 
Inelastic Solids by a New Stress Based Finite Element Method" 


p. 17, last line: 

replace "n (x, x)" by "n(X,T) M . 


p. 20, equation (3.26): 
replace " 0 " by 


p. 27, equation (3.52): 

replace -ji" by £ . 


P- 


43, 8 th line from page bottom: 

replace "In any case" by "In the case of solids without a 
natural time" 

5th line from page bottom 

replace "so the choice" by "so then the choice". 


p. 91, last line: 

replace "v" by "w" . 

p. 225, 3rd line from bottom of page: 

replace "p^" ^3 *x^" and "x^", respectively. 
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